Algorithm-Induced Prior for Image Restoration

I Introduction

Alternating direction methods of multipliers (ADMM) is perhaps the most popular algorithm for solving linear inverse problems in recent years, particularly for image restoration . Despite different perspectives of the algorithm, (e.g., operator splitting , proximal methods , split Bregman , to name a few,) the common principle behind ADMM is to convert the original minimization of the form

and solve the constrained problem by alternatingly minimizing the augmented Lagrangian function. Under mild conditions, e.g., when $f(\cdot)$ is strongly convex and $s(\cdot)$ is convex, the convergence of the algorithm is typically guaranteed .

In setting up the optimization problem in (1), the objective function $f(\cdot)$ and the regularization function $s(\cdot)$ are almost always fixed before running the algorithm. For example, when solving a non-blind deblurring problem using a total variation regularization , the functions $f(\cdot)$ and $s(\cdot)$ are

where $\boldsymbol{A}$ is the blur operator, $\boldsymbol{y}$ is the observed image, and $\|\boldsymbol{x}\|_{TV}$ is the total variation norm of the image $\boldsymbol{x}$.

For most image restoration problems, $f(\boldsymbol{x})$ is chosen according to the forward imaging model, and is fixed as long as we agree with the forward model. But $s(\boldsymbol{x})$ is the user’s subjective belief of how the solution should look like, a.k.a. the prior. In literature, apart from the total variation prior mentioned in (3), there are enormous number of priors we can use. However, there is one thing in common, which is that $s(\boldsymbol{x})$ has to be defined before using the ADMM algorithm.

In this paper, I present an ADMM algorithm where the regularization function $s(\boldsymbol{x})$ is unknown a-priori. At a first glance, this might seem unnatural because if $s(\boldsymbol{x})$ is unknown, then it is unclear about what we are trying to optimize in (1). However, as will be discussed shortly, the ADMM algorithm can generally be written as two modules – an inverse module, and a denoising module. The idea is to replace the denoising module by some off-the-shelf image denoising algorithm, e.g., non-local means or BM3D . In other words, we do not explicitly define $s(\boldsymbol{x})$ before running the algorithm, but use a denoising algorithm to perform the role of $s(\boldsymbol{x})$.

Replacing the denoising module of the ADMM algorithm by an off-the-shelf denoising algorithm was first proposed by Bouman and colleagues , to the best of my knowledge. Perhaps of the heuristic nature of the method, they call it the “Plug-and-Play” algorithm to stress that one can plug in any denoising algorithm and get the ADMM algorithm running. Under appropriate conditions on the denoising algorithm, one can prove the convergence of the plug-and-play .

In the context of compressive sensing , a similar version of the plug-and-play is also being studied. In , Baranuik and colleagues considered an approximate message passing (AMP) algorithm for recovering images. Recognizing that AMP also has a “inverse module” and a “denoising module”, they replace the shrinkage step in the conventional AMP with an off-the-shelf denoising algorithm (BM3D in their paper). Again, under appropriate conditions of the denoising algorithm, they proved that the AMP converges.

I-B Contributions

The focus of this paper is not to find weaker conditions under which “plug-and-play” converges. Rather, I like to address another equally important question: What is the original prior $s(\boldsymbol{x})$ if we choose a particular denoising algorithm? Answering this problem is essential to understand this type of algorithms in general. To make the discussion concrete, I will focus on the class of symmetric smoothing filters which is broad enough to include many denoising methods such as bilateral filter, non-local means and LARK , but at the same time also allows us to exploit matrix structures, e.g., the graph Laplacian . I call the new prior as an algorithm-induced prior to reflect the algorithmic nature of the prior.

The rest of the paper is organized as follows. I will first setup the problem in Section II. Then, in Section III, I will address the question about the original prior of the ADMM-induced algorithm, and discuss linkages with the conventional graph Laplacian prior. Experimental results are presented in Section IV, and a conclusion is given in Section V.

II Concept of Algorithm-Induced Prior

To begin the discussion I will first briefly introduce the ADMM algorithm. Interested readers can read for additional technical details.

Given the constrained minimization (2), the ADMM algorithm defines the augmented Lagranian function as

The minimizations in (5) and (6) are known as the primal updates, whereas the descent step in (7) is the dual update. If both $f(\cdot)$ and $s(\cdot)$ are closed, proper and convex, and if $\mathcal{L}(\cdot)$ has a saddle point, then one can prove convergence of the ADMM algorithm in terms of primal residue, primal objective and dual variable . In case when (5) and (6) are solved simultaneously instead of sequentially as presented above, then one will obtain the augmented Lagrangian method (ALM).

With some manipulations and rearrangement of terms we can show the following.

where $\boldsymbol{\bar{u}}^{(k)}\overset{\text{def}}{=}\frac{1}{\rho}\boldsymbol{u}^{(k)}$ is the scaled 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)}$.

The proof is skipped because it is essentially completing squares. The reason of rewriting the ADMM as above is to demonstrate the modular structure of the ADMM algorithm which we shall discuss shortly. As a side remark, the iterations (8)-(10) can be defined as a proximal algorithm .

To gain more insights into the modular structure presented in Proposition 1, let us consider the following example.

If we use $f(\cdot)$ and $s(\cdot)$ given in (3), we observe that (8) and (9) become

which is a reconstruction problem with a quadratic regularization, and a denoising problem (to denoise $\boldsymbol{\widetilde{v}}^{(k)}$) with a total variation regularization, respectively.

II-B Algorithm-Induced Prior

Recognizing the “denoising” module in the ADMM algorithm, we replace the $\boldsymbol{v}$-subproblem by a denoising algorithm. Formally, if we denote $\mathcal{D}_{h}$ as the denoising algorithm, i.e.,

Consider the non-local means as an example. One can first construct a kernel matrix $\boldsymbol{K}_{h}^{(k)}$ with the $(i,j)$-th entry

where $\boldsymbol{\widetilde{v}}^{(k)}_{i}$ denotes the $i$-th patch of the input $\boldsymbol{\widetilde{v}}^{(k)}$. Then, by applying Sinkhorn-Knopp balancing algorithm one can determine a pair of diagonal matrices $\boldsymbol{R}^{(k)}$ and $\boldsymbol{C}^{(k)}$ such that

III Analysis of Algorithm-Induced Prior

In this section I will address the question: What is the original prior $s(\boldsymbol{x})$ if we choose $\mathcal{D}_{h}$ as a symmetric smoothing filter? For notational simplicity I will drop the scripts $(\cdot)^{(k)}$ and $(\cdot)_{h}$.

The first main result is stated in the following proposition, which provides an explicit formula for the regularization $s(\boldsymbol{x})$ when symmetric smoothing filters are used.

There are two ways of proving this proposition. The first way is a “reverse engineering” approach. By plugging (14) into (15) and setting the first order derivative to zero we can show that $\boldsymbol{\widehat{v}}=\boldsymbol{W}\boldsymbol{\widetilde{v}}$.

The alternative proof is a constructive one. First, we observe that in order to obtain $\boldsymbol{W}\boldsymbol{\widetilde{v}}$ on the right hand side of (15), we must have $s(\boldsymbol{v})$ being quadratic. Therefore, we let

for some symmetric matrix $\boldsymbol{C}$ and constant $\alpha$. Taking the first order derivative of the resulting function yields

Rearranging the terms, we obtain a linear equation

Since $\alpha$ can be arbitrary, we set $\alpha=\rho/(2\lambda)$. Consequently, we have $(\boldsymbol{C}+\boldsymbol{I})\boldsymbol{v}=\boldsymbol{\widetilde{v}}$. Multiplying both sides with $\boldsymbol{W}$ yields $\boldsymbol{W}(\boldsymbol{C}+\boldsymbol{I})\boldsymbol{v}=\boldsymbol{W}\boldsymbol{\widetilde{v}}$. Thus, in order to obtain $\boldsymbol{v}=\boldsymbol{W}\boldsymbol{\widetilde{v}}$, $\boldsymbol{C}$ must be chosen such that

which gives $\boldsymbol{C}=(\boldsymbol{I}-\boldsymbol{W})\boldsymbol{W}^{+}$. ∎

There are a few important properties of the regularization $s(\boldsymbol{v})$ shown in (14). First, for any fixed $\boldsymbol{W}$, the matrix $(\boldsymbol{I}-\boldsymbol{W})\boldsymbol{W}^{+}$ is symmetric positive semidefinite. In fact, since the eigenvalues of a symmetric smoothing filter $\boldsymbol{W}$ is always bounded between 0 and 1, i.e., $0\preceq\boldsymbol{\Sigma}\preceq 1$, it holds that the eigenvalues of $(\boldsymbol{I}-\boldsymbol{W})\boldsymbol{W}^{+}$ is also bounded between 0 and 1. Therefore, if $\boldsymbol{W}$ is pre-defined before running the ADMM algorithm, then $s(\boldsymbol{x})$ is convex and hence the overall optimization is also convex.

Second, if we compare (14) with the conventional graph Laplacian regularization $s(\boldsymbol{v})=\boldsymbol{v}^{T}\boldsymbol{L}\boldsymbol{v}$ in the literature , where $\boldsymbol{L}\overset{\text{def}}{=}\boldsymbol{I}-\boldsymbol{W}$, we observe that (14) has an additional term $\boldsymbol{W}^{+}$. Using a graph signal processing terminology, we can view $\boldsymbol{W}$ as a lowpass filter and $\boldsymbol{L}$ is a highpass filter. $\boldsymbol{W}^{+}$ is a bandpass filter because of the truncation property of the pseudo-inverse. Therefore, the regularization $\boldsymbol{v}^{T}(\boldsymbol{I}-\boldsymbol{W})\boldsymbol{W}^{+}\boldsymbol{v}$ penalizes a smaller (but more focused) set of graph frequencies than the conventional regularization $\boldsymbol{v}^{T}\boldsymbol{L}\boldsymbol{v}$. In Section IV we will compare the performance.

III-B Closed-form Solution

Proposition 2 suggests a new prior which deserves a closer look. First of all, assume, for simplicity, that the matrix $\boldsymbol{W}$ is fixed throughout the ADMM iteration. This can be done either in an oracle setting (i.e., find $\boldsymbol{W}$ from the ground truth solution), or in a pre-filtering setting (i.e., find $\boldsymbol{W}$ from some initial guess of the solution). Both ways are common in image restoration .

Substituting $f(\boldsymbol{x})=\frac{1}{2}\|\boldsymbol{A}\boldsymbol{x}-\boldsymbol{y}\|^{2}$ and the specific prior $s(\boldsymbol{x})$ given by (14) into the original optimization (1), the problem becomes

which is a quadratic optimization. Closed-form solution of (16) exists, and is given by solving the normal equation

Since closed-form solution exists, it is possible to bypass the ADMM iterations and obtain the solution efficiently. However, from a computational perspective, there are two issues of (17) which we need to overcome. First, (17) involves inverting an $n\times n$ matrix which is computationally prohibitive for large $n$. Second, if $\boldsymbol{W}$ has a full rank but with some very small eigenvalues, $\boldsymbol{W}^{+}$ will cause numerical instability, depending on the numerical threshold for truncating the eigenvalues. Therefore, if we want to use the closed form solution in (17), one possible approach is to bypass the pseudo-inverse $\boldsymbol{W}^{+}$. This can be done using the following algebraic trick.

where $\boldsymbol{W}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{U}^{T}$ is the eigen-decomposition of $\boldsymbol{W}$.

Define $\boldsymbol{p}\overset{\text{def}}{=}\boldsymbol{W}^{+}\boldsymbol{x}$ (or, equivalently, $\boldsymbol{x}=\boldsymbol{W}\boldsymbol{p}$). Then it holds that

because $\boldsymbol{W}=\boldsymbol{W}^{T}$. Consider the eigen-decomposition $\boldsymbol{W}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{U}^{T}$, and let $\boldsymbol{q}=\boldsymbol{U}^{T}\boldsymbol{p}$, it follows that

The minimizer of this quadratic function is given by the solution of the normal equation

Since $\boldsymbol{U}\boldsymbol{U}^{T}=\boldsymbol{I}$, it holds that $\boldsymbol{p}=\boldsymbol{U}\boldsymbol{q}$ and hence the solution is $\boldsymbol{x}=\boldsymbol{W}\boldsymbol{p}=\boldsymbol{W}\boldsymbol{U}\boldsymbol{q}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{q}$. ∎

Inspecting (19), we observe that the matrix inversion only involves an $m\times m$ matrix, which is significantly smaller than the $n\times n$ matrix in (17). The matrix $\boldsymbol{A}\boldsymbol{V}$ are usually not difficult to evaluate. Below are two examples for image inpainting and deblurring.

For image inpainting, the matrix $\boldsymbol{A}$ is a binary diagonal matrix. Hence, $\boldsymbol{A}\boldsymbol{V}$ involves picking the non-zero columns of $\boldsymbol{V}$.

For image deblurring, the matrix $\boldsymbol{A}$ is a convolution. Hence, the multiplication of $\boldsymbol{A}$ and $\boldsymbol{v}_{i}$, the $i$th column of $\boldsymbol{V}$, is a blurring operation on $\boldsymbol{v}_{i}$.

In practice, symmetric smoothing filters sometimes have very narrow spatial support and hence $\boldsymbol{W}$ is a banded diagonal matrix. A banded diagonal $\boldsymbol{W}$ has a significantly higher rank, making the eigen-decomposition difficult. However, the good news is that such $\boldsymbol{W}$ is often easy to compute. In this case, the closed form should be replaced by the ADMM iteration.

IV Experimental Results

where $\boldsymbol{\eta}\sim\mathcal{N}(0,\sigma^{2})$ is an additive iid Gaussian noise. In this experiment, $\sigma=0.05$.

There are two choices of the filter $\boldsymbol{W}$. The first choice is to compute $\boldsymbol{W}$ from the ground truth solution. This is called the oracle setting, and is the best possible setting we can use under our framework. The second choice is to compute $\boldsymbol{W}$ from some initial estimate of the solution. In this problem, the initial estimate is performed using the classical Shepard’s interpolation method . If a more sophisticated initial estimator is used, it is likely that the performance will be improved.

We compare two graph Laplacian priors, namely

where $\boldsymbol{L}\overset{\text{def}}{=}\boldsymbol{I}-\boldsymbol{W}$ is the classical graph Laplacian, and $\boldsymbol{C}\overset{\text{def}}{=}(\boldsymbol{I}-\boldsymbol{W})\boldsymbol{W}^{+}$ is the proposed algorithm-induced prior. The results are shown in Table I. It is evident from the result that $s_{\boldsymbol{C}}(\boldsymbol{x})$ performs consistently better than $s_{\boldsymbol{L}}(\boldsymbol{x})$ for both the estimated $\boldsymbol{W}$ and the oracle $\boldsymbol{W}$. The gap is especially prominent if we look at the oracle case at 80% missing.

Figure 1 shows a visual comparison of an image captured by an i-Phone 6 camera with 50% missing pixel generated by MATLAB simulation. It should be reminded that in all experiments the parameter $\rho$ are adjusted accordingly for $s_{\boldsymbol{L}}(\boldsymbol{x})$ and $s_{\boldsymbol{C}}(\boldsymbol{x})$. Figure 1(d) illustrates such dependence: The optimal $\rho$ are different for different priors. However, the best PSNR of $s_{\boldsymbol{C}}(\boldsymbol{x})$ is significantly higher than that of $s_{\boldsymbol{L}}(\boldsymbol{x})$.

V Conclusion

Algorithm-induced prior is a strong performing but intriguing prior that we have little understanding about. Therefore, being able to explicitly write down the formula of the algorithm-induced prior is an important step which allows us to analyze the performance of such prior. In this paper, I demonstrated the case of symmetric smoothing filters and drew connections with the conventional graph Laplacian prior. On a set of image inpainting experiments, algorithm-induced prior offers consistently better results than the conventional graph Laplacian. As we progress along this direction, I believe that the interplay between the objective function and the denoising procedure should be studied in greater details.

VI Acknowledgement

I like to thank Charles Bouman for pointing me to the problem, and Suhas Sreehari for fruitful discussions. I also thank Dror Baron for telling me his work on AMP using BM3D.

References