Plug-and-Play ADMM for Image Restoration: Fixed Point Convergence and Applications

I Introduction

for some conditional probability $p(\boldsymbol{y}\,|\,\boldsymbol{x})$ defining the forward imaging model, and a prior distribution $p(\boldsymbol{x})$ defining the probability distribution of the latent image. Because of the explicit use of the forward and the prior models, MAP estimation is also a model-based image reconstruction (MBIR) method [Bouman_2015] which has many important applications in deblurring [Afonso_Bioucas-Dias_Figueiredo_2010, Chan_Khoshabeh_Gibson_2011, Yang_Zhang_Yin_2009], interpolation [Dahl_Hansen_Jensen_2010, Garcia_2010, Zhou_Chen_Ren_2009], super-resolution [Dong_Loy_He_2014, Peleg_Elad_2014, He_Siu_2011, Yang_Wright_Huang_2008] and computed tomography [Sreehari_Venkatakrishnan_Wohlberg_2015], to name a few.

It is not difficult to see that solving the MAP problem in (1) is equivalent to solving an optimization problem

with $f(\boldsymbol{x})\overset{\text{def}}{=}-\log p(\boldsymbol{y}\,|\,\boldsymbol{x})$ and $g(\boldsymbol{x})\overset{\text{def}}{=}-(1/\lambda)\log p(\boldsymbol{x})$. The optimization in (2) is a generic unconstrained optimization. Thus, standard optimization algorithms can be used to solve the problem. In this paper, we focus on the alternating direction method of multiplier (ADMM) [Boyd_Parikh_Chu_Peleato_Eckstein_2011], which has become the workhorse for a variety of problems in the form of (2).

The idea of ADMM is to convert (2), an unconstrained optimization, into a constrained problem

and consider its augmented Lagrangian function:

The minimizer of (3) is then the saddle point of $\mathcal{L}$, which can be found by solving a sequence of subproblems

where $\boldsymbol{\bar{u}}^{(k)}\overset{\text{def}}{=}(1/\rho)\boldsymbol{u}^{(k)}$ is the scaled Lagrange multiplier, $\boldsymbol{\widetilde{x}}^{(k)}\overset{\text{def}}{=}\boldsymbol{v}^{(k)}-\boldsymbol{\bar{u}}^{(k)}$ and $\boldsymbol{\widetilde{v}}^{(k)}\overset{\text{def}}{=}\boldsymbol{x}^{(k+1)}+\boldsymbol{\bar{u}}^{(k)}$. Under mild conditions, e.g., when both $f$ and $g$ are closed, proper and convex, and if a saddle point of $\mathcal{L}$ exists, one can show that the iterates (5)-(7) converge to the solution of (3) (See [Boyd_Parikh_Chu_Peleato_Eckstein_2011] for details).

I-B Plug-and-Play ADMM

An important feature of the ADMM iterations (5)-(7) is its modular structure. In particular, (5) can be regarded as an inversion step as it involves the forward imaging model $f(\boldsymbol{x})$, whereas (6) can be regarded as a denoising step as it involves the prior $g(\boldsymbol{v})$. To see the latter, if we define $\sigma=\sqrt{\lambda/\rho}$, it is not difficult to show that (6) is

Treating $\boldsymbol{\widetilde{v}}^{(k)}$ as the “noisy” image, (8) minimizes the residue between $\boldsymbol{\widetilde{v}}^{(k)}$ and the “clean” image $\boldsymbol{v}$ using the prior $g(\boldsymbol{v})$. For example, if $g(\boldsymbol{x})=\|\boldsymbol{x}\|_{TV}$ (the total variation norm), then (8) is the standard total variation denoising problem.

Building upon this intuition, Venkatakrishnan et al. [Venkatakrishnan_Bouman_Wohlberg_2013] proposed a variant of the ADMM algorithm by suggesting that one does not need to specify $g$ before running the ADMM. Instead, they replace (6) by using an off-the-shelf image denoising algorithm, denoted by $\mathcal{D}_{\sigma}$, to yield

Because of the heuristic nature of the method, they called the resulting algorithm as the Plug-and-Play ADMM. An interesting observation they found in [Venkatakrishnan_Bouman_Wohlberg_2013] is that although Plug-and-Play ADMM appears ad-hoc, for a number of image reconstruction problems the algorithm indeed performs better than some state-of-the-art methods. A few recent reports have concurred similar observations [Sreehari_Venkatakrishnan_Wohlberg_2015, Dar_Bruckstein_Elad_2015, Rond_Giryes_Elad_2015, Brifman_Romano_Elad_2016].

I-C Challenges of Plug-and-Play ADMM

From a theoretical point of view, the main challenge of analyzing Plug-and-Play ADMM is the denoiser $\mathcal{D}_{\sigma}$. Since $\mathcal{D}_{\sigma}$ is often nonlinear and does not have closed form expressions, the analysis has been very difficult. Specifically, the following three questions remain open:

Convergence of the Algorithm. Classical results of ADMM require $g$ to be closed, proper and convex in order to ensure convergence [Boyd_Parikh_Chu_Peleato_Eckstein_2011]. While newer results have extended ADMM for nonconvex problems [Hong_Luo_Razaviyayn_2015], there is little work addressing the case when $g$ is defined implicitly through $\mathcal{D}_{\sigma}$. To the best of our knowledge, the only existing convergence analysis, to date, is the one by Sreehari et al. [Sreehari_Venkatakrishnan_Wohlberg_2015] for the case when $\mathcal{D}_{\sigma}$ is a symmetric smoothing filter [Milanfar_2013b, Chan_Zickler_Lu_2015]. However, for general $\mathcal{D}_{\sigma}$ the convergence is not known.

Original Prior. Since $\mathcal{D}_{\sigma}$ is an off-the-shelf image denoising algorithm, it is unclear what prior $g$ does it correspond to. In [Chan_2016], Chan addresses this question by explicitly deriving the original prior $g$ when $\mathcal{D}_{\sigma}$ is a symmetric smoothing filter [Chan_2016]. In this case, the author shows that $g$ is a modified graph Laplacian prior, with better restoration performance compared to the conventional graph Laplacian [Milanfar_2013a]. However, beyond symmetric smoothing filters it becomes unclear if we can find the corresponding $g$.

Implementation. The usage of Plug-and-Play ADMM has been reported in a few scattered occasions, with some work in electron tomography [Sreehari_Venkatakrishnan_Wohlberg_2015], compressive sensing [Dar_Bruckstein_Elad_2015], and some very recent applications in Poisson recovery [Rond_Giryes_Elad_2015] and super-resolution [Brifman_Romano_Elad_2016]. However, the common challenge underpinning these applications is whether one can obtain a fast solver for the inversion step in (5). This has not been a problem for conventional ADMM, because in many cases we can use another variable splitting strategy to replace $\boldsymbol{v}=\boldsymbol{x}$ in (3), e.g., using $\boldsymbol{v}=\boldsymbol{B}\boldsymbol{x}$ when $g(\boldsymbol{x})=\|\boldsymbol{B}\boldsymbol{x}\|_{1}$ [Afonso_Bioucas-Dias_Figueiredo_2010].

I-D Related Works

Plug-and-Play ADMM was first reported in 2013. Around the same period of time there is an independent series of studies using denoisers for approximate message passing (AMP) algorithms [Metzler_Maleki_Baraniuk_2014, Metzler_Maleki_Baraniuk_2015, Tan_Ma_Baron_2014, Tan_Ma_Baron_2015]. The idea was to replace the shrinkage step of the standard AMP algorithm with any off-the-shelf algorithm in the class of “proper denoisers” – denoisers which ensure that the noise variance is sufficiently suppressed. (See Section II-C for more discussions.) However, this type of denoise-AMP algorithms rely heavily on the Gaussian statistics of the random measurement matrix $\boldsymbol{A}$ in a specific forward model $f(\boldsymbol{x})=\|\boldsymbol{A}\boldsymbol{x}-\boldsymbol{y}\|^{2}$. Thus, if $f(\boldsymbol{x})$ departs from quadratic or if $\boldsymbol{A}$ is not random, then the behavior of the denoise-AMP becomes unclear.

Using denoisers as building blocks of an image restoration algorithm can be traced back further, e.g., wavelet denoiser for signal deconvolution [Neelamani_Choi_Baraniuk_2004]. Of particular relevance to Plug-and-Play ADMM is the work of Danielyan et al. [Danielyan_Katkovnik_Egiazarian_2012], where they proposed a variational method for deblurring using BM3D as a prior. The idea was later extended by Zhang et al. to other restoration problems [Zhang_Zhao_Gao_2014]. However, these algorithms are customized for the specific denoiser BM3D. In contrast, the proposed Plug-and-Play ADMM supports any denoiser satisfying appropriate assumptions. Another difference is that when BM3D is used in [Danielyan_Katkovnik_Egiazarian_2012] and [Zhang_Zhao_Gao_2014], the grouping of the image patches are fixed throughout the iterations. Plug-and-Play ADMM allows re-calculation of the grouping at every iteration. In this aspect, the Plug-and-Play ADMM is more general than these algorithms.

A large number of denoisers we use nowadays are patch-based denoising algorithms. All these methods can be considered as variations in the class of universal denoisers [Weissman_Ordentlich_Seroussi_2005, Sivaramakrishnan_Weissman_2009] which are asymptotically optimal and do not assume external knowledge of the latent image (e.g., prior distribution). Asymptotic optimality of patch-based denoisers has been recognized empirically by Levin et al. [Levin_Nadler_2011, Levin_Nadler_Durand_2012], who showed that non-local means [Buades_Coll_2005_Journal] approaches the MMSE estimate as the number of patches grows to infinity. Recently, Ma et al. [Ma_Zhu_Baron_2016] made attempts to integrate universal denoisers with approximate message passing algorithms.

I-E Contributions

The objective of this paper is to address the first and the third issue mentioned in Section I-C. The contributions of this paper are as follows:

First, we modify the original Plug-and-Play ADMM by incorporating a continuation scheme. We show that the new algorithm is guaranteed to converge for a broader class of denoisers known as the bounded denoisers. Bounded denoisers are asymptotically invariant in the sense that the denoiser approaches an identity operator as the denoising parameter vanishes. Bounded denoisers are weaker than the non-expansive denoisers presented in [Sreehari_Venkatakrishnan_Wohlberg_2015]. However, for weaker denoisers we should also expect a weaker form of convergence. We prove that the new Plug-and-Play ADMM has a fixed point convergence, which complements the global convergence results presented in [Sreehari_Venkatakrishnan_Wohlberg_2015].

Second, we discuss fast implementation techniques for image super-resolution and single photon imaging problems. For the super-resolution problem, conventional ADMM requires multiple variable splits or an inner conjugate gradient solver to solve the subproblem. We propose a polyphase decomposition based method which gives us closed-form solutions. For the single photon imaging problem, existing ADMM algorithm are limited to explicit priors such as total variation. We demonstrate how Plug-and-Play ADMM can be used and we present a fast implementation by exploiting the separable feature of the problem.

The rest of the paper is organized as follows. We first discuss the Plug-and-Play ADMM algorithm and the convergence properties in Section II. We then discuss the applications in Section III. Experimental results are presented in Section IV.

II Plug-and-Play ADMM and Convergence

In this section we present the proposed Plug-and-Play ADMM and discuss its convergence property. Throughout this paper, we assume that the unknown image $\boldsymbol{x}$ is bounded in an interval $[x_{\min},\;x_{\max}]$ where the upper and lower limits can be obtained from experiment or from prior knowledge. Thus, without loss of generality we assume $\boldsymbol{x}\in^{n}$.

The proposed Plug-and-Play ADMM algorithm is a modification of the conventional ADMM algorithm in (5)-(7). Instead of choosing a constant $\rho$, we increase $\rho$ by $\rho_{k+1}=\gamma_{k}\rho_{k}$ for $\gamma_{k}\geq 1$. In optimization literature, this is known as a continuation scheme [Ng_2002] and has been used in various problems, e.g., [Harmany_Marcia_Willet_2011, Wang_Yang_Yin_2008]. Incorporating this idea into the ADMM algorithm, we obtain the following iteration:

where $\mathcal{D}_{\sigma_{k}}$ is a denoising algorithm (called a “denoiser” for short), and $\sigma_{k}\overset{\text{def}}{=}\sqrt{\lambda/\rho_{k}}$ is a parameter controlling the strength of the denoiser.

There are different options in setting the update rule for $\rho_{k}$. In this paper we present two options. The first one is a monotone update rule which defines

for a constant $\gamma>1$. The second option is an adaptive update rule by considering the relative residue:

For any $\eta\in[0,\,1)$ and let $\gamma>1$ be a constant, we conditionally update $\rho_{k}$ according to the followings:

If $\Delta_{k+1}\geq\eta\Delta_{k}$, then $\rho_{k+1}=\gamma\rho_{k}$.

If $\Delta_{k+1}<\eta\Delta_{k}$, then $\rho_{k+1}=\rho_{k}$.

The adaptive update scheme is inspired from [Goldstein_ODonoghue_Setzer_2012], which was originally used to accelerate ADMM algorithms for convex problems. It is different from the residual balancing technique commonly used in ADMM, e.g., [Boyd_Parikh_Chu_Peleato_Eckstein_2011], as $\Delta_{k+1}$ sums of all primal and dual residues instead of treating them individually. Our experience shows that the proposed scheme is more robust than residual balancing because the denoiser could potentially generate nonlinear effects to the residuals. Algorithm 1 shows the overall Plug-and-Play ADMM.

In the original Plug-and-Play ADMM by Venkatakrishnan et al. [Venkatakrishnan_Bouman_Wohlberg_2013], the update scheme is $\rho_{k}=\rho$ for some constant $\rho$. This is valid when the denoiser $\mathcal{D}_{\sigma}$ is non-expansive and has symmetric gradient. However, for general denoisers which could be expansive, the update scheme for $\rho_{k}$ becomes crucial to the convergence. (See discussion about non-expansiveness in Section II-B.)

Many denoising algorithms nowadays such as BM3D and non-local means require one major parameter A denoising algorithm often involves many other “internal” parameters. However, as these internal parameters do not have direct interaction with the ADMM algorithm, in this paper we keep all internal parameters in their default settings to simplify the analysis., typically an estimate of the noise level, to control the strength of the denoiser. In our algorithm, the parameter $\sigma_{k}$ in (11) is reminiscent to the noise level. However, unlike BM3D and non-local means where $\sigma_{k}$ is directly linked to the standard deviation of the i.i.d. Gaussian noise, in Plug-and-Play ADMM we treat $\sigma_{k}$ simply as a tunable knob to control the amount of denoising because the residue $(\boldsymbol{v}-\boldsymbol{\widetilde{v}}^{(k)})$ at the $k$th iterate is not exactly Gaussian. The adoption of the Gaussian denoiser $\mathcal{D}_{\sigma_{k}}$ is purely based on the formal equivalence between (8) and a Gaussian denoising problem.

In this paper, we assume that the parameter $\lambda$ is pre-defined by the user and is fixed. Its role is similar to the regularization parameter in the conventional ADMM problem. Tuning $\lambda$ can be done using external tools such as cross validation [Nguyen_Milanfar_Golub_2001] or SURE [Ramani_Blu_Unser_2008].

II-B Global and Fixed Point Convergence

Before we discuss the convergence behavior, we clarify two types of convergence.

We refer to the type of convergence in the conventional ADMM as global convergence, i.e., convergence in primal residue, primal objective and dual variables. To ensure global convergence, one sufficient condition is that $g$ is convex, proper and closed [Boyd_Parikh_Chu_Peleato_Eckstein_2011]. For Plug-and-Play ADMM, a sufficient condition is that $\mathcal{D}_{\sigma}$ has symmetric gradient and is non-expansive [Sreehari_Venkatakrishnan_Wohlberg_2015]. In this case, $g$ exists due to a proximal mapping theorem of Moreau [Moreau_1965]. However, proving non-expansive denoisers could be difficult as it requires

for any $\boldsymbol{x}$ and $\boldsymbol{y}$, with $\kappa\leq 1$. Even for algorithms as simple as non-local means, one can verify numerically that there exists pairs $(\boldsymbol{x},\boldsymbol{y})$ that would cause $\kappa>1$. In the Appendix we demonstrate a counter example.

Since $\mathcal{D}_{\sigma}$ can be arbitrary and we do not even know the existence of $g$, we consider fixed point convergence instead. Fixed point convergence guarantees that a nonlinear algorithm can enter into a steady state asymptotically. In nonlinear dynamical systems, these limit points are referred to as the stable-fixed-points. For any initial guess lying in a region called the basin of attraction the algorithm will converge [Wiggins_1990]. For Plug-and-Play ADMM, we conjecture that fixed point convergence is the best we can ask for unless further assumptions are made on the denoisers.

II-C Convergence Analysis of Plug-and-Play ADMM

We define the class of bounded denoisers.

for some universal constant $C$ independent of $n$ and $\sigma$.

Bounded denoisers are asymptotically invariant in the sense that it ensures $\mathcal{D}_{\sigma}\rightarrow\mathcal{I}$ (i.e., the identity operator) as $\sigma\rightarrow 0$. It is a weak condition which we expect most denoisers to have. The asymptotic invariant property of a bounded denoiser prevents trivial mappings from being considered, e.g., $\mathcal{D}_{\sigma}(\boldsymbol{x})=0$ for all $\boldsymbol{x}$.

It would be useful to compare a bounded denoiser with a “proper denoiser” defined in [Metzler_Maleki_Baraniuk_2014]. A proper denoiser $\widetilde{\mathcal{D}}_{\sigma}$ is a mapping that denoises a noisy input $\boldsymbol{x}+\sigma\boldsymbol{\epsilon}$ with the property that

for any $\kappa<1$, where $\boldsymbol{\epsilon}\sim\mathcal{N}(0,\boldsymbol{I})$ is the i.i.d. Gaussian noise. Note that in (17), we require the input to be a deterministic signal $\boldsymbol{x}$ plus an i.i.d. Gaussian noise. Moreover, the parameter must match with the noise level. In contrast, a bounded denoiser can take any input and any parameter.

Besides the conditions on $\mathcal{D}_{\sigma}$ we also assume that the negative log-likelihood function $f$ has bounded gradients:

The main convergence result of this paper is as follows.

(Fixed Point Convergence of Plug-and-Play ADMM). Under Assumption 1 and for any bounded denoiser $\mathcal{D}_{\sigma}$, the iterates of the Plug-and-Play ADMM defined in Algorithm 1 demonstrates a fixed-point convergence. That is, there exists $(\boldsymbol{x}^{*},\boldsymbol{v}^{*},\boldsymbol{u}^{*})$ such that $\|\boldsymbol{x}^{(k)}-\boldsymbol{x}^{*}\|_{2}\rightarrow 0$, $\|\boldsymbol{v}^{(k)}-\boldsymbol{v}^{*}\|_{2}\rightarrow 0$ and $\|\boldsymbol{u}^{(k)}-\boldsymbol{u}^{*}\|_{2}\rightarrow 0$ as $k\rightarrow\infty$.

Intuitively, what Theorem 1 states is that as $k\rightarrow\infty$, the continuation scheme forces $\rho_{k}\rightarrow\infty$. Therefore, the inversion in (10) and the denoising in (11) have reducing influence as $\rho_{k}$ grows. Hence, the algorithm converges to a fixed point. Theorem 1 also ensures that $\boldsymbol{x}^{(k)}\rightarrow\boldsymbol{v}^{(k)}$ which is an important property of the original Plug-and-Play ADMM algorithm [Sreehari_Venkatakrishnan_Wohlberg_2015]. The convergence of $\boldsymbol{x}^{(k)}\rightarrow\boldsymbol{v}^{(k)}$ holds because $\boldsymbol{u}^{(k+1)}=\boldsymbol{u}^{(k)}+(\boldsymbol{x}^{(k+1)}-\boldsymbol{v}^{(k+1)})$ converges. In practice, experimentally we observe that if the algorithm is terminated early to reduce the runtime, then $\boldsymbol{v}^{(k)}$ tends to provide a slightly better solution.

II-D Stopping Criteria

Since we are seeking for fixed point convergence, a natural stopping criteria is to determine if $\|\boldsymbol{x}^{(k+1)}-\boldsymbol{x}^{(k)}\|_{2}$, $\|\boldsymbol{v}^{(k+1)}-\boldsymbol{v}^{(k)}\|_{2}$ and $\|\boldsymbol{u}^{(k+1)}-\boldsymbol{u}^{(k)}\|_{2}$ are sufficiently small. Following the definition of $\Delta_{k+1}$ in (15), we choose to terminate the iteration when

for some tolerance level $\mathtt{tol}$. Alternatively, we can also terminate the algorithm when

where $\epsilon_{1}=\|\boldsymbol{x}^{(k+1)}-\boldsymbol{x}^{(k)}\|_{2}/\sqrt{n}$, $\epsilon_{2}=\|\boldsymbol{v}^{(k+1)}-\boldsymbol{v}^{(k)}\|_{2}/\sqrt{n}$ and $\epsilon_{3}=\|\boldsymbol{u}^{(k+1)}-\boldsymbol{u}^{(k)}\|_{2}/\sqrt{n}$.

In practice, the tolerance level does not need to be extremely small in order to achieve good reconstruction quality. In fact, for many images we have tested, setting $\mathtt{tol}\approx 10^{-3}$ is often sufficient. Figure 1 provides a justification. In this experiment, we tested an image super-resolution problem for 10 testing images (See Configuration 3 in Section IV-A for details). It can be observed that the PSNR becomes steady when $\mathtt{tol}$ drops below $10^{-3}$. Moreover, size of the image does not seem to be an influencing factor. Smaller images such as Cameraman256, House256 and Peppers256 shows similar characteristics as bigger images. The more influencing factor is the combination of the update ratio $\gamma$ and the initial value $\rho_{0}$. However, unless $\gamma$ is close to 1 and $\rho_{0}$ is extremely small (which does not yield good reconstruction anyway), our experience is that setting $\mathtt{tol}$ at $10^{-3}$ is usually valid for $\gamma\in(1,2)$ and $\rho_{0}\in(10^{-5},10^{-2})$.

The choice of the initial parameter $\rho_{0}$ requires some tuning but is typically good for $\rho_{0}\in(10^{-5},10^{-2})$. Figure 2 shows the behavior of the algorithm for different values of $\rho_{0}$, ranging from $10^{0}$ to $10^{-4}$. We compare the original Plug-and-Play ADMM (i.e., with constant $\rho_{k}=\rho_{0}$, the red lines), monotone update rule (i.e., $\rho_{k+1}=\gamma\rho_{k}$, the blue lines), and the adaptive update rule (the black lines). We make two observations regarding the difference between the proposed algorithm and the original Plug-and-Play ADMM [Sreehari_Venkatakrishnan_Wohlberg_2015]:

Stability: The original Plug-and-Play ADMM [Sreehari_Venkatakrishnan_Wohlberg_2015] requires a highly precise $\rho_{0}$. For example, in Figure 2 the best PSNR is achieved when $\rho_{0}=1$; When $\rho_{0}$ is less than $10^{-2}$, the PSNR becomes very poor. The proposed algorithm works for a much wider range of $\rho_{0}$.

Final PSNR: The proposed Plug-and-Play ADMM is a generalization of the original Plug-and-Play ADMM. The added degrees of freedom are the new parameters $(\rho_{0},\gamma,\eta)$. The original Plug-and-Play ADMM is a special case when $\gamma=1$. Therefore, for optimally tuned parameters, the proposed Plug-and-Play ADMM is always better than or equal to the original Plug-and-Play ADMM. This is verified in Figure 2, which shows that the best PSNR is attained by the proposed method.

II-F Initial Guesses

The initial guesses $\boldsymbol{x}^{(0)}$, $\boldsymbol{v}^{(0)}$ and $\boldsymbol{u}^{(0)}$ have less impact to the final PSNR. This can be seen from Figure 3. In this experiment, we randomly draw 100 initial guesses $\boldsymbol{x}^{(0)}$ from a uniform distribution in $^{n}$. The auxiliary variable is set as $\boldsymbol{v}^{(0)}=\boldsymbol{x}^{(0)}$, and the Lagrange multiplier $\boldsymbol{u}^{(0)}$ is 0. As shown in Figure 3, the initial guesses do not cause significant difference in term of PSNR at the limit. The standard deviation at the limit is 0.0059 dB, implying that with 99.7% probability (3 standard deviations) the PSNR will stay within $\pm 0.0176$ dB from its average.

III Applications

As we discussed in the introduction, Plug-and-Play ADMM algorithm has a wide range of applications. However, in order to enable the denoising step, Plug-and-Play ADMM uses a specific variable splitting strategy. The challenge it brings, therefore, is whether we can solve the subsequent subproblems efficiently. The purpose of this section is to address this issue by presenting two applications where fast implementation can be achieved.

Image super-resolution can be described by a linear forward model with two operations: an anti-aliasing filter and a subsampling process. The function $f(\boldsymbol{x})$ is quadratic in the form

Consequently, the solution is the pseudo-inverse

For special cases of $\boldsymbol{H}$ and $\boldsymbol{S}$, (21) has known efficient implementation as follows.

Non-blind deblurring is a special case when $\boldsymbol{S}=\boldsymbol{I}$. In this case, since $\boldsymbol{H}$ is circulant which is diagonalizable by the discrete Fourier transform matrices, (21) can be efficiently implemented by

where $\mathcal{F}(\cdot)$ is the Fourier transform operator, $h$ is the finite impulse response filter representing the blur kernel, $\overline{(\cdot)}$ is the complex conjugate, and the multiplication/division are element-wise operations.

Image interpolation is a special case when $\boldsymbol{H}=\boldsymbol{I}$. In this case, since $\boldsymbol{S}^{T}\boldsymbol{S}$ is a diagonal matrix with binary entries, (21) can be efficiently implemented using an element-wise division:

III-B Polyphase Implementation for Image Super-Resolution

When $\boldsymbol{G}=\boldsymbol{S}\boldsymbol{H}$, solving the $f$-subproblem becomes non-trivial because $\boldsymbol{H}^{T}\boldsymbol{S}^{T}\boldsymbol{S}\boldsymbol{H}$ is neither diagonal nor diagonalizable by the Fourier transform. In literature, the two most common approaches are to introduce multi-variable split to bypass (21) (e.g., [Afonso_Bioucas-Dias_Figueiredo_2010, Almeida_Figueiredo_2013]) or use an inner conjugate gradient to solve (21) (e.g., [Brifman_Romano_Elad_2016]). However, multi-variable splitting requires additional Lagrange multipliers and internal parameters. It also generally leads to slower convergence than single-variable split. Inner conjugate gradient is computationally expensive as it requires an iterative solver. In what follows, we show that when $\boldsymbol{S}$ is the standard $K$-fold downsampler (i.e., sub-sample the spatial grid uniformly with a factor $K$ along horizontal and vertical directions), and when $\boldsymbol{H}$ is a circular convolution, it is possible to derive a closed-form solution We assume the boundaries are circularly padded. In case of other types boundary conditions or unknown boundary conditions, we can pre-process the image by padding the boundaries circularly. Then, after the super-resolution algorithm we crop the center region. The alternative approach is to consider multiple variable split as discussed in [Almeida_Figueiredo_2013]. .

Our closed form solution begins by considering the Sherman-Morrison-Woodbury identity, which allows us to rewrite (21) as

The more critical step is the following observation. We note that the matrix $\boldsymbol{G}\boldsymbol{G}^{T}$ is given by

Since $\boldsymbol{S}$ is a $K$-fold downsampling operator, $\boldsymbol{S}^{T}$ is a $K$-fold upsampling operator. Defining $\boldsymbol{\widetilde{H}}=\boldsymbol{H}\boldsymbol{H}^{T}$, which can be implemented as a convolution between the blur kernel $h$ and its time-reversal, we observe that $\boldsymbol{S}\boldsymbol{\widetilde{H}}\boldsymbol{S}^{T}$ is a “upsample-filter-downsample” sequence. This idea is illustrated in Figure 4.

We next study the polyphase decomposition [Vaidyanathan_1992] of Figure 4. Polyphase decomposition allows us to write

where $\widetilde{H}(z)$ is the $z$-transform representation of the blur matrix $\boldsymbol{\widetilde{H}}=\boldsymbol{H}\boldsymbol{H}^{T}$, and $\widetilde{H}_{k}(z^{K})$ is the $k$th polyphase component of $\widetilde{H}(z)$. Illustrating (25) using a block diagram, we show in Figure 5 the decomposed structure of Figure 4. Then, using Noble identity [Vaidyanathan_1992], the block diagram on the left hand side of Figure 5 becomes the one shown on the right hand side. Since for any $k>1$, placing a delay $z^{-k}$ between an upsampling and a downsampling operator leads to a zero, the overall system simplifies to a finite impulse response filter $\widetilde{H}_{0}(z)$, which can be pre-computed.

We summarize this by the following proposition.

The operation of $\boldsymbol{S}\boldsymbol{H}\boldsymbol{H}^{T}\boldsymbol{S}^{T}$ is equivalent to applying a finite impulse response filter $\widetilde{H}_{0}(z)$, which is the 0th polyphase component of the filter $\boldsymbol{H}\boldsymbol{H}^{T}$.

To implement the 0th polyphase component, we observe that it can be done by downsampling the convolved filter $\boldsymbol{\widetilde{H}}=\boldsymbol{H}\boldsymbol{H}^{T}$. This leads to the procedure illustrated in Algorithm 2.

The implication of Proposition 1 is that since $\boldsymbol{G}\boldsymbol{G}^{T}$ is equivalent to a finite impulse response filter $\widetilde{h}_{0}$, (24) can be implemented in closed-form using the Fourier transform:

where we recall that $\boldsymbol{b}=\boldsymbol{G}^{T}\boldsymbol{y}+\rho\boldsymbol{\widetilde{x}}$.

The effectiveness of the proposed closed-form solution can be seen from Figure 6. In this figure, we compare with a brute force conjugate gradient method presented in [Brifman_Romano_Elad_2016]. When $\boldsymbol{H}$ satisfies periodic boundary conditions, the closed-form solution is exact. If the boundaries are not periodic, alternative solutions can be considered, e.g., [Matakos_Ramani_Fessler_2013].

III-C Application 2: Single Photon Imaging

The second application is a single photon imaging problem using quanta image sensors (QIS) [Fossum_2011]. Using ADMM for QIS was previously reported in [Chan_Lu_2014, Elgendy_Chan_2016]. Here, we show how the Plug-and-Play ADMM can be used for the problem.

and $\alpha$ is a sensor gain. Given $\boldsymbol{s}$, the photons arriving at the sensors follow a Poisson distribution with a rate given by $\boldsymbol{s}$. Let $Z_{i}$ be the random variable denoting the number of photons at jot $i$, we have

The final QIS output, $Y_{i}$, is a binary bit resulted from truncating $Z_{i}$ using a threshold $q$. That is,

When $q=1$, the probability of observing $Y_{i}=y_{i}$ given $s_{i}$ is

The recovery goal is to estimate $\boldsymbol{x}$ from the observed binary bits $\boldsymbol{y}$. Taking the negative log and summing over all pixels, the function $f$ is defined as

where $K_{j}^{1}=\sum_{i=1}^{K}y_{(j-1)K+i}$ is the number of ones in the $j$th unit pixel, and $K_{j}^{0}=\sum_{i=1}^{K}(1-y_{(j-1)K+i})$ is the number of zeros in the $j$th unit pixel. (Note that for any $j$, $K^{1}_{j}+K^{0}_{j}=K$.) Consequently, substituting (29) into (10) yields the $f$-subproblem

Since this optimization is separable, we can solve each individual variable $x_{j}$ independently. Thus, for every $j$, we solve a single-variable optimization by taking derivative with respect to $x_{j}$ and setting to zero, yielding

which is a one-dimensional root finding problem. By constructing an offline lookup table in terms of $K_{0}$, $\rho$ and $\widetilde{x}_{j}$, we can solve (30) efficiently.

IV Experimental Results

In this section we present the experimental results. For consistency we use BM3D in all experiments, although other bounded denoisers will also work. We shall not compare Plug-and-Play ADMM using different denoisers as it is not the focus of the paper.

We consider a set of 10 standard test images for this experiment as shown in Figure 7. All images are gray-scaled, with sizes between $256\times 256$ and $512\times 512$. Four sets of experimental configurations are studied, and are shown in Table I.