Solving RED with Weighted Proximal Methods

TL;DR

This paper introduces weighted proximal methods (WPMs) to improve the computational efficiency of the RED framework by reducing the number of denoiser calls needed for inverse problem solving.

Contribution

It demonstrates that existing RED solvers are special cases of WPMs and proposes variants that significantly decrease computation time.

Findings

WPMs encompass existing RED solvers as special cases.

Variants of WPM reduce denoiser calls and computation time.

Numerical experiments confirm improved efficiency of WPM-based RED solutions.

Abstract

REgularization by Denoising (RED) is an attractive framework for solving inverse problems by incorporating state-of-the-art denoising algorithms as the priors. A drawback of this approach is the high computational complexity of denoisers, which dominate the computation time. In this paper, we apply a general framework called weighted proximal methods (WPMs) to solve RED efficiently. We first show that two recently introduced RED solvers (using the fixed point and accelerated proximal gradient methods) are particular cases of WPMs. Then we show by numerical experiments that slightly more sophisticated variants of WPM can lead to reduced run times for RED by requiring a significantly smaller number of calls to the denoiser.

Tables

Table 1. TABLE I: Denoiser evaluations required to attain a similar PSNR. The first and second rows per each image refer to image deblurring and the third row refers to super-resolution. The minimal number of denoiser evaluations is marked in bold.

	FP-MPE	APG	WPM
Butterfly	$54$	$34$	$\mathrm{𝟐𝟓}$
	$54$	$26$	$\mathrm{𝟏𝟕}$
	$80$	$51$	$\mathrm{𝟐𝟔}$
Boats	$24$	$\mathrm{𝟐𝟎}$	$21$
	$60$	$34$	$\mathrm{𝟐𝟐}$
	$36$	$20$	$\mathrm{𝟏𝟐}$
House	$24$	$\mathrm{𝟏𝟖}$	$19$
	$62$	$26$	$\mathrm{𝟐𝟓}$
	$18$	$15$	$\mathrm{𝟏𝟎}$
Parrot	$39$	$30$	$\mathrm{𝟐𝟎}$
	$52$	$40$	$\mathrm{𝟑𝟔}$
	$49$	$31$	$\mathrm{𝟐𝟖}$
Lena	$48$	$34$	$\mathrm{𝟐𝟗}$
	$47$	$16$	$\mathrm{𝟏𝟓}$
	$37$	$26$	$\mathrm{𝟏𝟖}$
Barbara	$14$	$12$	$\mathrm{𝟏𝟏}$
	$48$	$23$	$\mathrm{𝟏𝟔}$
	$17$	$15$	$\mathrm{𝟏𝟏}$
Peppers	$42$	$29$	$\mathrm{𝟐𝟐}$
	$41$	$40$	$\mathrm{𝟑𝟒}$
	$38$	$30$	$\mathrm{𝟐𝟖}$
Leaves	$50$	$41$	$\mathrm{𝟑𝟒}$
	$36$	$18$	$\mathrm{𝟏𝟒}$
	$60$	$41$	$\mathrm{𝟏𝟐}$

Equations

$${\mathbf{x}}_{ML}^{\ast }=\arg \underset{\mathbf{x}}{\min }\frac{1}{2}\Vert \mathbf{H}\mathbf{x}-\mathbf{y}{\Vert }_2^2.$$

$$R(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\mathcal{T}}\left(\mathbf{x}-f(\mathbf{x})\right),$$

$$\begin{array}[]{rcl}{\bm{x}}^{*}_{MAP}&=&\arg\max_{\bm{x}}P(\bm{x}|\bm{y})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(MAP)\\ &=&\arg\max_{\bm{x}}P(\bm{y}|\bm{x})P(\bm{x})~{}~{}~{}~{}~{}~{}~{}~{}(Bayes~{}rule)\\ &=&\arg\min_{\bm{x}}-\log\{P(\bm{y}|\bm{x})\}-\log{P(\bm{x})}.\end{array}$$

$$P(\mathbf{x})=constant\cdot e^{-\alpha R(\mathbf{x})},$$

$$P(\mathbf{y}|\mathbf{x})=constant\cdot e^{-\frac{1}{2{\sigma }^2}\Vert \mathbf{H}\mathbf{x}-\mathbf{y}{\Vert }_2^2},$$

$${\mathbf{x}}_{MAP}^{\ast }=\arg \underset{\mathbf{x}}{\min }E(\mathbf{x})\triangleq \frac{1}{2{\sigma }^2}\Vert \mathbf{H}\mathbf{x}-\mathbf{y}{\Vert }_2^2+\alpha R(\mathbf{x}).$$

$${\mathbf{x}}_{MAP}^{\ast }=\arg \underset{\mathbf{x}}{\min }E(\mathbf{x})\triangleq \frac{1}{2{\sigma }^2}\Vert \mathbf{H}\mathbf{x}-\mathbf{y}{\Vert }_2^2+\frac{\alpha }{2}{\mathbf{x}}^{\mathcal{T}}\left(\mathbf{x}-f(\mathbf{x})\right).$$

$${\nabla }_{x}f(\mathbf{x})\mathbf{x}=\underset{ϵ\to 0}{\lim }\frac{f(\mathbf{x}+ϵ\mathbf{x})-f(\mathbf{x})}{ϵ}=f(\mathbf{x}).$$

$${\nabla }_{\mathbf{x}}R(\mathbf{x})=\mathbf{x}-f(\mathbf{x}).$$

$${\nabla }_{\mathbf{x}}E(\mathbf{x})=\frac{1}{{\sigma }^2}{\mathbf{H}}^{\mathcal{T}}\left(\mathbf{H}\mathbf{x}-\mathbf{y}\right)+\alpha \left(\mathbf{x}-f(\mathbf{x})\right).$$

$$\frac{1}{{\sigma }^2}{\mathbf{H}}^{\mathcal{T}}\left(\mathbf{H}{\mathbf{x}}_{k+1}-\mathbf{y}\right)+\alpha \left({\mathbf{x}}_{k+1}-f({\mathbf{x}}_k)\right)=0.$$

$$\underset{\mathbf{x}}{\min }\mathcal{F}(\mathbf{x})\triangleq g(\mathbf{x})+h(\mathbf{x}),$$

$${\text{prox}}_{h,\mathbf{B}}(\hat{\mathbf{x}})=\arg \underset{\mathbf{u}}{\min }\left\{h(\mathbf{u})+\frac{1}{2}\Vert \mathbf{u}-\hat{\mathbf{x}}{\Vert }_{\mathbf{B}}^2\right\},$$

$${\nabla }_{\mathbf{u}}h(\mathbf{u})+\mathbf{B}(\mathbf{u}-\hat{\mathbf{x}})=0.$$

$$\left(\frac{a_k}{{\sigma }^2}{\mathbf{H}}^{\mathcal{T}}\mathbf{H}+{\mathbf{B}}_k\right){\mathbf{x}}_{k+1}=\frac{a_k}{{\sigma }^2}{\mathbf{H}}^{\mathcal{T}}\mathbf{y}+{\mathbf{B}}_k{\mathbf{x}}_k-a_k\alpha \left({\mathbf{x}}_k-f({\mathbf{x}}_k)\right).$$

Full text

I Introduction

THE goal of inverse problems is to recover an unknown signal $\bm{x}\in\Re^{N}$ from an indirect measurement $\bm{y}\in\Re^{M}$. The measurement is commonly modelled as $\bm{y}={\mathcal{H}}(\bm{x})+\bm{e},$ where $\mathcal{H}(\cdot)$ denotes an abstract operator and $\bm{e}$ is often assumed to be white Gaussian noise with mean zero and variance $\sigma^{2}$. In this paper, we assume $\mathcal{H}(\cdot)$ to be a linear operator, $\mathcal{H}(\bm{x})=\bm{H}\bm{x},$ with $\bm{H}\in\Re^{M\times N}$, and focus on natural images. Lacking any prior knowledge about the signal $\bm{x}$, we may reconstruct $\bm{x}$ via the maximum likelihood (ML) minimization problem,

[TABLE]

However, it is well-known that this approach is not generally useful. Even in the simple denoising problem, where $\bm{H}$ is the identity matrix, ML results in $\bm{x}^{*}_{ML}=\bm{y}$, that is, we simply recover the noisy image. Furthermore, quite often $M<N$, resulting in infinitely many solutions, and even if this is not the case, $\bm{H}$ may be highly ill-conditioned. For these reasons, the prevalent approach is to assume that the signal $\bm{x}$ is sampled from some prior distribution, and to employ the maximum a posteriori probability (MAP) estimator, as formulated in Section II. In our setting, MAP will result in adding to the right-hand side of (1) a term $\alpha R(\bm{x})$, where $\alpha$ is a parameter and $R$ is a regularization as discussed below. This approach has been applied with a large variety of priors, such as $\ell_{2}$-based regularization [1], wavelets [2], total variation [3], kernel regularization [4], sparsity [5], and neural networks [6].

Naturally, the most widely studied problem in this framework is image denoising, e.g., [7, 5, 8, 6, 9]. Indeed, recent work suggests that the performance of leading image denoisers is close to a possible ceiling [10, 11, 12]. The availability of such powerful denoising algorithms has motivated researchers to seek ways to employ denoisers as priors for quite general inverse problems. The authors in [13, 14, 15] “manually” adopted priors used in existing denoisers for specific alternative inverse problems. Following this path, several authors proposed a general framework, called Plug-and-Play Priors ($P^{3}$) [16, 17], for using the abundance of high-performance image denoisers as priors for other inverse problems. These authors formulate inverse problems as an optimization task and employ an Alternating Direction Method of Multipliers (ADMM) algorithm to tackle the corresponding minimization problem [18]. The image denoising algorithm is incorporated in each step of ADMM as an implicit prior.

Motivated by $P^{3}$, Romano et al. introduced REgularization by Denoising (RED) [19], which defines an optimization problem that includes the denoiser as an explicit prior. Given a differentiable denoiser $f(\bm{x})$, RED employs the following prior,

[TABLE]

where $\cdot^{\mathcal{T}}$ denotes transpose. Under two assumptions this leads to a convex minimization problem, and standard gradient based iterative methods are guaranteed to converge to a global minimum. Further details are provided in Section II.

Using state-of-the-art denoisers to construct priors is appealing, as it enables us to exploit the vast progress in denoising algorithms for addressing general inverse problems, and RED is a good framework to achieve this goal due to its flexibility. However, RED may be relatively expensive because at each iteration we must apply the denoising algorithm to evaluate the gradient, and the complexity of denoising algorithms is generally high. Indeed, the numerical experiments in [19] reveal this concern. In that paper the authors propose three solvers for RED, namely, steepest descent (SD), the fixed-point (FP) method and the ADMM scheme. Amongst these, the FP method is the most efficient, but it still needs hundreds of iterations to complete the recovery process.

Recently, [20] employed vector extrapolation to accelerate the FP method for RED, whereas [21] applies an accelerated proximal gradient (APG) algorithm111APG is also known as FISTA [22] or Nesterov’s acceleration [23].. Both these approaches are faster than FP for RED, but they still require dozens of iterations. In this paper, we propose a general framework called weighted proximal methods (WPMs) [24]. We show that FP and APG are in fact two particular variants of WPMs, and that by seeking a more effective weighting for WPMs we obtain a faster algorithm.

The rest of this paper is organized as follows. We review the RED framework and the FP and APG solvers in Section II. The general WPM scheme is introduced in Section III, and the choice of weighting is discussed. Numerical experiments on image deblurring and super-resolution tasks are presented in Section IV to demonstrate the efficiency of WPMs, followed by conclusions in Section V.

II REgularization by Denoising (RED)

The MAP recovery process is formulated as follows:

[TABLE]

Assuming a robust Gibbs-like distribution of $\bm{x}$, we have

[TABLE]

where $R$ denotes the so-called prior and $\alpha>0$ is a scaling parameter. Note that small $R(\bm{x})$ corresponds to highly probable signals. If $\bm{e}$ is sampled from white Gaussian noise with mean zero and variance $\sigma^{2}$, then we have

[TABLE]

This leads to the following minimization problem [25],

[TABLE]

Substituting the RED prior (2) into (3), we obtain

[TABLE]

In [19] two assumptions are made regarding the image denoising algorithm used in RED:

Assumption 1**.**

For any scalar $c$ arbitrarily close to $1$, $f(c\bm{x})=cf(\bm{x})$.

Assumption 2**.**

The spectral radius of the symmetric Jacobian $\nabla_{\bm{x}}f(\bm{x})$ is upper bounded by one.

Given Assumption 1, we have

[TABLE]

Hence, the gradient of $R(\bm{x})$ is the residual of the denoiser,

[TABLE]

With (5), the gradient of $E(\bm{x})$ becomes

[TABLE]

Assumption 2 implies convexity of $R(\bm{x})$, and therefore of $E(\bm{x})$ as well. Hence, any solution of $\nabla_{\bm{x}}E(\bm{x})=0$ yields a global minimum. This is a nonlinear problem, and we therefore resort to iterative solvers. One such solver is the FP method mentioned above, which lags the nonlinear term $f(\bm{x})$:

[TABLE]

We note that (7) can efficiently be solved for $\bm{x}_{k+1}$ exactly in the Fourier domain if $\bm{H}$ is block-circulant, or treated iteratively for a general $\bm{H}$. The FP method can be accelerated using the APG approach as described in the following algorithm. Further discussion of APG can be found in [21].

III Weighted Proximal Methods

Consider the following composite problem and assume its solution set is nonempty,

[TABLE]

where $g$ and $h$ are convex and differentiable. Denote the proximal operator by

[TABLE]

where $\bm{B}$ is a symmetric positive definite matrix called the weighting and $\|\cdot\|_{\bm{B}}$ denotes the $\bm{B}$-norm, $\|\bm{q}\|_{\bm{B}}=\sqrt{\bm{q}^{\mathcal{T}}\bm{B}\bm{q}}$. With these, we describe the explicit form of WPMs for (8) in Algorithm 2 [24, Chap. 10.7.5]. Note that by setting $\bm{B}=\beta\bm{I}$ with $\beta>0$, we recover the proximal gradient (PG) method. Usually, PG is used for (8) when $h$ is nonsmooth [26], whereas here we use it even though $h$ is differentiable. We do this for computational efficiency, knowing that applying the denoiser is the most expensive part of the solution process.

To apply Algorithm 2 to RED, we set $g(\bm{x})=\alpha R(\bm{x})=\frac{\alpha}{2}\bm{x}^{\mathcal{T}}\left(\bm{x}-f(\bm{x})\right)$ and $h(\bm{x})=\frac{1}{2\sigma^{2}}\|\bm{H}\bm{x}-\bm{y}\|_{2}^{2}$. If $h(\bm{x})$ is convex, solving (9) is equivalent to satisfying the first-order optimality condition,

[TABLE]

Substituting $\hat{\bm{x}}\leftarrow\bm{x}_{k}-a_{k}\bm{B}_{k}^{-1}\nabla_{\bm{x}}g(\bm{x}_{k})$, $h(\bm{u})\leftarrow a_{k}h(\bm{u})$ and $\bm{u}\leftarrow\bm{x}_{k+1}$, at the $k$th iteration into (10) and rearranging, we obtain

[TABLE]

In this paper, we use the conjugate gradient (CG) method to approximately solve (11) for $\bm{x}_{k+1}$.

Next we discuss possible practical choices for the weighting $\bm{B}_{k}$. Note first that if we set $\bm{B}_{k}=\alpha\bm{I}$, where $\bm{I}$ is the identity matrix, and select the step-size $a_{k}=1$, (11) is reduced to (7) and we recover the FP method. Moreover, by using the accelerated version of Algorithm 2 (cf. [24, Chap. 10.7.5]) we get APG [21] . We now propose a more elaborate approach of choosing some approximation to the Hessian of $\alpha R(\bm{x})$ as the weighting. (Because of the abstract denoiser in $R(\bm{x})$, the exact Hessian is not computable.) Specifically, we choose the symmetric-rank-one (SR1) approximation to the Hessian [27, Chap. 6.2], as is used in quasi-Newton methods. The SR1 approximation is described in Algorithm 3. This choice yields faster convergence in our experiments than either FP or APG, as shown below. We henceforth use WPM to denote Algorithm 2 with the weighting chosen by Algorithm 3.

Unlike the traditional SR1, we formulate each $\bm{B}_{k}$ from the initial $\bm{H}_{0}$ rather the previous iterate $\bm{B}_{k-1}$ [27]. Moreover, we scale $\bm{H}_{0}$ by $\gamma>1$ as suggested in [28], which we found useful in practice. In the practical implementation of Algorithm 3, we efficiently represent $\bm{B}_{k}$ as a matrix-vector multiplication operator rather than as an explicit matrix.

In general, the step-size $a_{k}$ in Algorithm 2 needs to be chosen by some line search process to guarantee monotonically decreasing objective values at each iteration. However, because evaluating the objective value in RED requires calling the denoiser, standard line search methods may dramatically increase the complexity of the algorithm. To maintain a low computational cost, we fix $a_{k}=1$ and reduce the step-size by half only if the objective value exhibits a relative growth above some threshold, i.e., $E(\bm{x}_{k+1})-E(\bm{x}_{k})>\varepsilon E(\bm{x}_{k+1})$, where we use $\varepsilon=10^{-2}$ in all our experiments. In practice, we found that we never needed to reduce the step-size.

In this paper we only investigate the SR1 approximation to the Hessian of $\alpha R(\bm{x})$. We acknowledge that a more accurate Hessian estimate may prove to be even more cost-effective for RED, but leave such investigation to future work. Because we use an approximate Hessian for the weighting, our algorithm is equivalent to a quasi-newton proximal method. It follows that if both $g(\bm{x})$ and $h(\bm{x})$ are strongly convex and their gradients are Lipschitz continuous, WPM with SR1 estimation, an appropriate step-size $a_{k}$, and exact solution of (9), converges linearly; see details in [29]. Because we depart from these strict requirements for efficiency, we cannot claim provable convergence in our implementation. However, in all our experiments WPM converged. Finally, we note that [21] challenges the validity in practice of the underlying assumptions of RED for most denoisers, concluding that (6) is not truly the gradient of (4). Nevertheless, setting (6) to zero, as is the objective of all the algorithms we discuss here, remains a most attractive method for signal recovery.

IV Numerical Experiments

In this section we investigate the performance of solvers for RED. Following [19], we perform our tests on image deblurring and super-resolution tasks and use the trainable nonlinear reaction diffusion (TNRD) [6] method as the abstract denoiser. We remark that one can adopt deep denoising techniques instead of TNRD, since the differentiability requirement of the denoiser is not mandatory in practice [21]. This may possibly lead to improved results in practice, but we do not investigate such options here. Also, since the authors in [19] already show the superiority of RED for image deblurring and super-resolution tasks compared with other popular algorithms, we largely omit such comparisons in this paper and concentrate on computational efficiency. Moreover, the experiments conducted in [20] demonstrated that the FP method converges faster than LBFGS and Nesterov’s acceleration for RED. Therefore, we only compare WPM to FP [19], FP-MPE [20], and APG [21]. All of the experiments are carried out on a laptop with Intel i$7-6500$U CPU @2.50GHz and 8GB RAM.

For image deblurring, the image is degraded by convolving with a point spread function (PSF), $9\times 9$ uniform blur or a Gaussian blur with a standard derivation $1.6$, and then adding Gaussian noise with mean zero and $\sigma=\sqrt{2}$. The recovered peak-signal-to-noise ratio (PSNR) versus the number of denoiser evaluations (left column) and running time (right column) when using RED for the “Starfish” image are shown in Figure 2. We find that the performances of FP-MPE and APG are similar, whereas WPM is more efficient than both, requiring less denoiser evaluations and running time to achieve a comparable PSNR. These results also indicate that indeed the denoiser dominates the complexity of solving RED.

Next, we test the algorithms on image super-resolution. A low resolution image is generated by blurring a high-resolution image with a $7\times 7$ Gaussian kernel with standard derivation $1.6$, and then downscaling by a factor of 3. To the resulting image we add Gaussian noise with mean zero and $\sigma=5$, resulting in our deteriorated image. The PSNR of the recovered fine-resolution image versus the number of denoiser evaluations (left) and running time (right) for the “Plants” image are presented in Figure 3. Again, we observe that WPM requires less denoiser evaluations and running time to achieve a comparable PSNR.

Examining the performance of the algorithms further, we run them on eight additional images tested in [19]. For each image, we run the FP method with $200$ denoiser evaluations and take the final PSNR as a benchmark. Then we examine how many denoiser evaluations are needed for APG, FP-MPE and WPM, to achieve a similar PSNR. The results are listed in Table I. Evidently, with the exception of “Boats” and “House” in the deblurring task, we observe that WPM requires the smallest number of denoiser evaluations to achieve a comparable PSNR, demonstrating its efficiency for solving RED. Additionally, we present the recovered results of the “Starfish” and “Leaves” images from deblurring with uniform and Gaussian blurs, respectively, and the “Butterfly” image from super-resolution in Figure 1 to visualize the effectiveness of RED solved by WPM.

V Conclusion

In this paper, we propose a general framework for RED called weighted proximal methods (WPMs). By setting $\bm{B}_{k}=\alpha\bm{I}$ and $a_{k}=1$, we retrieve the FP and APG methods. However, by choosing the weighting to be an approximation to the Hessian of $\alpha R(\bm{x})$, we obtain a more efficient algorithm. The experiments on image deblurring and super-resolution tasks demonstrate that WPM with a simple and inexpensive approximation to the Hessian can substantially reduce the overall number of denoiser evaluations in the recovery process, usually resulting in significant speedup. In future work we aim to design better Hessian approximations in order to accelerate the computation further.