Synthesis versus analysis in patch-based image priors

TL;DR

This paper explores the analysis and synthesis approaches in patch-based image priors, establishing their differences, proposing a synthesis model, and demonstrating its practical application in image denoising.

Contribution

It introduces a synthesis formulation for patch-based image priors, contrasting it with existing analysis models, and provides an ADMM-based method for denoising.

Findings

All existing patch-based models are analysis-based.

A new synthesis formulation is proposed and formalized.

The synthesis approach is computationally feasible for denoising.

Abstract

In global models/priors (for example, using wavelet frames), there is a well known analysis vs synthesis dichotomy in the way signal/image priors are formulated. In patch-based image models/priors, this dichotomy is also present in the choice of how each patch is modeled. This paper shows that there is another analysis vs synthesis dichotomy, in terms of how the whole image is related to the patches, and that all existing patch-based formulations that provide a global image prior belong to the analysis category. We then propose a synthesis formulation, where the image is explicitly modeled as being synthesized by additively combining a collection of independent patches. We formally establish that these analysis and synthesis formulations are not equivalent in general and that both formulations are compatible with analysis and synthesis formulations at the patch level. Finally, we present an instance of the alternating direction method of multipliers (ADMM) that can be used to perform image denoising under the proposed synthesis formulation, showing its computational feasibility. Rather than showing the superiority of the synthesis or analysis formulations, the contributions of this paper is to establish the existence of both alternatives, thus closing the corresponding gap in the field of patch-based image processing.

Equations

$$\mathbf{y}=\mathbf{x}+\mathbf{w},$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP}\in {\text{mbox}argmax}_{\mathbf{x}}p(\mathbf{x}|\mathbf{y})={\text{mbox}argmin}_{\mathbf{x}}\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{x}{\Vert }_2^2-\log p_X(\mathbf{x}),$$

$$-\log p_X(\mathbf{x})=\phi (\mathbf{F}\mathbf{x})+b,$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP-A}\in {\text{mbox}argmax}_{\mathbf{x}}p(\mathbf{x}|\mathbf{y})={\text{mbox}argmin}_{\mathbf{x}}\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{x}{\Vert }_2^2+\phi (\mathbf{F}\mathbf{x}),$$

$$\hat{\boldsymbol{x}}_{\mbox{\scriptsize MAP-S}}=\boldsymbol{G}\hat{\mbox{\boldmath$\theta$}},\;\mbox{where}\;\hat{\mbox{\boldmath$\theta$}}\in\mathop{\mbox{argmin}}_{\mbox{\boldmath$\theta$}}\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{G}\mbox{\boldmath$\theta$}\|_{2}^{2}+\zeta(\mbox{\boldmath$\theta$}),$$

$$p_X(\mathbf{x})=\frac{1}{Z}\prod _{m=1}^{M}f({\mathbf{P}}_m\mathbf{x}),$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP-A}\in \arg \underset{\mathbf{x}}{\min }\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{x}{\Vert }_2^2-\sum _{m=1}^M\log f({\mathbf{P}}_m\mathbf{x}).$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP-A}$$

$${\hat{x}}_{\text{mbox}MAP-A}$$

$$\mathbf{P}\mathbf{x}=\left[\begin{array}{c}{\mathbf{P}}_1\mathbf{x}\\ ⋮\\ {\mathbf{P}}_M\mathbf{x}\end{array}\right]\in {\mathbb{R}}^{Mn}.$$

$$p_X(\mathbf{x})\propto p_V(\mathbf{P}\mathbf{x}),$$

$$p_V(\mathbf{v})=\prod _{m=1}^{M}f({\mathbf{v}}_m).$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP-A}\in \arg \underset{\mathbf{x}}{\min }\frac{1}{2{\sigma }^2}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2^2-\log p_V(\mathbf{P}\mathbf{x}),$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP-A}$$

$$\mathbf{x}=\sum _{m=1}^M{\mathbf{Q}}_m{\mathbf{z}}_m,$$

$${\mathbf{Q}}_1=\left[\begin{array}{cc}1/2& 0\\ 0& 1/2\\ 0& 0\\ 0& 0\end{array}\right],$$

$${\mathbf{Q}}_3=\left[\begin{array}{cc}0& 0\\ 0& 0\\ 1/2& 0\\ 0& 1/2\end{array}\right],$$

$$\mathbf{x}=\mathbf{Q}\mathbf{z}.$$

$$\log p_Z(\mathbf{z})=\sum _{m=1}^M\log g({\mathbf{z}}_m).$$

$${\hat{\mathbf{x}}}_{\text{mbox}MAP-S}=\mathbf{Q}\hat{\mathbf{z}},\text{mbox}where\hat{\mathbf{z}}\in \arg \underset{\mathbf{z}}{\min }\frac{1}{2{\sigma }^2}\Vert \mathbf{Q}\mathbf{z}-\mathbf{y}{\Vert }_2^2-\log p_Z(\mathbf{z}).$$

$$\hat{\boldsymbol{x}}_{\mbox{\scriptsize MAP-S}}=\boldsymbol{Q}\mbox{\boldmath$\Delta$}\hat{\mbox{\boldmath$\gamma$}},$$

$$\hat{\mbox{\boldmath$\gamma$}}\in\arg\min_{\mbox{\boldmath$\gamma$}}\frac{1}{2\sigma^{2}}\|\boldsymbol{Q}\mbox{\boldmath$\Delta$}\mbox{\boldmath$\gamma$}-\boldsymbol{y}\|_{2}^{2}-\log q_{\Gamma}(\mbox{\boldmath$\gamma$}).$$

$$\mathbf{Q}\mathbf{P}=\mathbf{I},$$

$$\mathbf{P}\mathbf{Q}\ne \mathbf{I}.$$

$$\begin{array}{rl}{\hat{\mathbf{x}}}_{\text{mbox}MAP-A}=\mathbf{Q}\hat{\mathbf{z}},\text{mbox}where\hat{\mathbf{z}}\in & \arg \underset{\mathbf{z}}{\min }\frac{1}{2{\sigma }^2}\Vert \mathbf{Q}\mathbf{z}-\mathbf{y}{\Vert }_2^2-\log p_V(\mathbf{z}),\\ & \text{mbox}subjectto\mathbf{z}\in \text{mbox}range(\mathbf{P}),\end{array}$$

$$\begin{array}{rl}& \underset{\mathbf{z},\mathbf{u}}{\min }\frac{1}{2{\sigma }^2}\Vert \mathbf{Q}\mathbf{z}-\mathbf{y}{\Vert }_2^2+\xi (\mathbf{u}),\\ & \text{mbox}subjectto\mathbf{z}=\mathbf{u},\end{array}$$

$$\xi (\mathbf{u})=-\log p_Z(\mathbf{u})=\sum _{m=1}^M\underset{{\xi }_m({\mathbf{u}}_m)}{\underbrace{(-\log g({\mathbf{u}}_m))}}$$

$${\mathbf{z}}^{(t+1)}$$

$${\mathbf{u}}^{(t+1)}$$

$${\mathbf{d}}^{(t+1)}$$

$${\mathbf{u}}_m^{(t+1)}={\text{mbox}argmin}_{{\mathbf{u}}_m}{\xi }_m({\mathbf{u}}_m)+\frac{\rho }{2}\Vert {\mathbf{z}}_m^{(t+1)}-{\mathbf{u}}_m-{\mathbf{d}}_m^{(t)}{\Vert }_2^2,$$

$$\boldsymbol{z}^{(t+1)}=\bigl{(}\boldsymbol{Q}^{T}\boldsymbol{Q}+\sigma^{2}\rho\boldsymbol{I}\bigr{)}^{-1}\bigl{(}\boldsymbol{Q}^{T}\boldsymbol{y}+\sigma^{2}\rho\,(\boldsymbol{u}^{(t)}+\boldsymbol{d}^{(t)})\bigr{)}.$$

Full text

1 Introduction

The use of patches in image processing can be seen as an instance of the “divide and conquer” principle: since it is admittedly very difficult to formulate a global prior/model for images, patch-based approaches use priors/models for patches (rather than whole images), the combination of which yields the desired image prior/model. To keep the discussion and formulation at their essential and focus on the image modelling aspects, we will concentrate on image denoising, arguably the quintessential image processing problem. Nevertheless, much of what will be presented below can be easily extended (at least in principle) to more general inverse problems.

There are basically two approaches to patch-based image denoising. In earlier methods [3], [5], [6], patches are extracted from the noisy image, then processed/denoised independently (or maybe even collaboratively, as in BM3D [5]), and finally returned to their original locations. Since the patches overlap (to avoid blocking artifacts), there are several estimates of each pixel, which are combined by some form of averaging (unweighted or weighted). This approach is also used in the nonlocal Bayesian method [9], and in methods based on Gaussian mixtures [13]. For a comprehensive review of these and related methods, see [9]. Arguably, a conceptual flaw of these methods is that they obtain patch estimates without explicitly taking into account that these will subsequently be combined additively. As a consequence, although some of these methods achieve state-of-the-art results, they do not explicitly provide a global image prior/model.

A more recent class of approaches does build global image models that are based on a function computed from image patches, but does not treat them as independent by explicitly taking into account that they are overlapping patches of the same image; this approach was initiated with the expected patch log-likelihood (EPLL) [15], and is adopted by most of the recent work [4], [12]. These methods do not have the conceptual flaw pointed out in the previous paragraph and provide a coherent global image model.

The analysis vs synthesis dichotomy in global image models/priors (e.g., based on wavelet frames, or total variation) has been first formalized in [7], and further studied in [11]; more recently, it has been ported to patch-wise models [4]. To the best of our knowledge, this dichotomy has not been pointed out before concerning the way in which patch-level models/priors are used to build a global image models; that is precisely the central contribution of this paper.

In the synthesis vs analysis dichotomy, the EPLL-type class of patch-based models can be seen as an analysis method (as explained below in detail). This paper shows that there exists the synthesis counterpart of EPLL; in other words, that the synthesis vs analysis dichotomy is also present in the way the whole image and the patches are related. We also show that the two formulations are not, in general, equivalent.

The remaining sections of this paper are organized as follows. After reviewing the classical analysis/synthesis dichotomy in Section 2, we shown in Section 3 that the classical patch-based methods follow an analysis formulation. Section 4 then introduces a synthesis patch-based formulation, and its relationship with the analysis counterpart is established in Section 5. In Section 6, we present an ADMM algorithm to efficiently perform image denoising under the proposed synthesis formulation. Finally, Section 7 concludes the paper by referring to future work directions.

2 Analysis vs Synthesis in Image Denoising

Before addressing patch-based models, we briefly review the analysis and synthesis global formulations of image denoising [7], where the goal is to estimate an unknown image $\boldsymbol{x}\in\mathbb{R}^{N}$ ($N$ is the total number of pixels in $\boldsymbol{x}$, which is a vectorized version of the corresponding $\sqrt{N}\times\sqrt{N}$ image) from a noisy version thereof

[TABLE]

where $\boldsymbol{w}$ is a sample of a white Gaussian noise field of zero mean and (known) variance $\sigma^{2}$, that is, $p(\boldsymbol{y}|\boldsymbol{x})=\mathcal{N}(\boldsymbol{y};\boldsymbol{x},\sigma^{2}\boldsymbol{I})$.

The classical approach to estimate $\boldsymbol{x}$ from $\boldsymbol{y}$ is to adopt (or learn) a prior $p_{X}$ for the unknown $\boldsymbol{x}$ and seek a maximizer of the posterior density (a maximum a posteriori–MAP–estimate)

[TABLE]

where $\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{x}\|_{2}^{2}=-\log p(\boldsymbol{y}|\boldsymbol{x})$, up to an irrelevant constant. The analysis and synthesis formulations build priors for $\boldsymbol{x}$ as follows.

Analysis:

in this formulation, the prior $p_{X}$ takes the form

[TABLE]

where $b$ is a constant and $\boldsymbol{F}$ a linear (analysis) operator. In this case, the MAP estimate takes the form

[TABLE]

where MAP-A stands for MAP-analysis.

Synthesis:

here, the starting point is to assume that $\boldsymbol{x}$ is linearly synthesized/represented according to \boldsymbol{x}=\boldsymbol{G}\mbox{\boldmath\theta}, where $\theta$ is a vector of coefficients, and the prior is formulated on $\theta$ rather than directly on $\boldsymbol{x}$; once an estimate \hat{}\mbox{\boldmath\theta} is obtained, the corresponding estimate of $\boldsymbol{x}$ is simply \hat{\boldsymbol{x}}=\boldsymbol{G}\hat{\mbox{\boldmath\theta}}. In summary,

[TABLE]

and \zeta(\mbox{\boldmath\theta})=-\log_{\Theta}(\mbox{\boldmath\theta}) is a negative log-prior on $\theta$.

One of the main distinguishing features of analysis and synthesis formulations is that in the former the object of the estimation procedure is the image itself, whereas in the latter, one estimates a representation from which the image estimate is synthesized.

3 Patch-Analysis Formulation

3.1 Expected Patch Log-likelihood

A central tool in most patch-based approaches is a collection of operators $\{\boldsymbol{P}_{m},\,m=1,...,M\}$ that extracts $M$ patches from a given image of size $\sqrt{N}\times\sqrt{N}$; each $\boldsymbol{P}_{m}$ can be seen as a binary matrix with size $n\times N$, where $n$ is the total number of pixels in each patch (assumed square, of size $\sqrt{n}\times\sqrt{n}$). The standard way of formulating a patch-based prior is by writing

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}_{+}$ is a function expressing the patch-wise prior distribution and $Z$ is a normalizing constant. Function $f$ may itself be a probability density function (pdf), e.g., a GMM, in which case this prior is an instance of a so-called product of experts (PoE [8]). However, $f$ does not need to be a pdf, as long as it takes non-negative values; in fact, (6) can also be seen as a factor graph model, where each factor corresponds to a patch and the all factors have the same function $f$ [14]. Moreover, this prior is also equivalent to the formulation known as EPLL (expected patch log-likleihood [15]), although EPLL was not originally interpreted as a prior.

Given a noisy image $\boldsymbol{y}$, a MAP estimate of $\boldsymbol{x}$ is given by

[TABLE]

3.2 Half Quadratic Splitting and ADMM

To tackle the large-scale optimization problem in (7), the so-called half quadratic splitting strategy replaces it with

[TABLE]

where $\boldsymbol{v}=(\boldsymbol{v}_{m},\;m=1,...,M)$, which obviously becomes equivalent to (7) as $\beta\rightarrow\infty$ [15]. The optimization problem in (8) is tackled by alternating between minimizing with respect to $\boldsymbol{v}$ and $\boldsymbol{x}$, while slowly increasing $\beta$. Other strategies for setting $\beta$ have also been proposed [15].

An obvious alternative (as recently mentioned in [10]) is to reformulate (7) as a constrained problem

[TABLE]

and tackle it with ADMM (alternating direction method of multipliers) [2]. Of course, convergence of ADMM for this problem can only be guaranteed if the negative log factors $-\log f$ are convex; this is not the case if $f$ is a GMM, but it is true if $-\log f$ is an $\ell_{r}$ norm or the $r$-th power thereof, e.g., $-\log f(\boldsymbol{v}_{m})=\|\boldsymbol{v}_{m}\|_{r}^{r}$, with $r\geq 1$.

Examining (8) (with $\beta\rightarrow\infty$) or (3.2) reveals that this formulation seeks a consensus among the patches, in the sense that the several replicates of each pixels that exist in different patches are forced to agree on a common value for that pixel. In other words, the clean patches are not modelled as additively generating a clean image, and are merely used to write a joint prior $p_{X}$ that factorizes across overlapping patches.

3.3 Identification as Analysis Formulation

The estimation criterion in (7) clearly falls in the analysis-type category [7], since it considers the image itself as the object to be estimated and as the argument of the prior. However, as a generative model for clean images, its meaning is not very clear, since it is not a trivial task to obtain samples from this distribution.

To obtain a more compact notation, let $\boldsymbol{P}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{Mn}$ be the operator (an $Mn\times N$ matrix) that extracts the set of $M$ patches, i.e.,

[TABLE]

The prior $p_{X}$ may then be written as

[TABLE]

where $p_{V}$ denotes a density defined in $\mathbb{R}^{Mn}$ according to

[TABLE]

With this notation, the MAP denoising problem (7) can be written as

[TABLE]

which clearly reveals the analysis nature of this formulation (see (4)).

3.4 Patch level Models

We stress that the analysis/synthesis dichotomy addressed in this paper concerns the way in which an image relates to its patches, not the way the patches themselves are modelled. In fact, the patch-analysis formulation just reviewed is compatible with a synthesis patch model, e.g., one that models each patch $\boldsymbol{P}_{m}\boldsymbol{x}$ as linear combination of elements of some dictionary $\boldsymbol{D}$, with coefficients \mbox{\boldmath\gamma}_{m} equipped with some prior $p_{\Gamma}$ (for example, a sparsity-inducing prior, as in [12]). Using the half quadratic splitting approach, the formulation becomes

[TABLE]

where \mbox{\boldmath\gamma}=(\mbox{\boldmath\gamma}_{m},\;m=1,...,M). Naturally, the patch-analysis formulation is also compatible with an analysis patch model [4], by using a patch prior of the form $f=\phi\circ\boldsymbol{A}$, i.e., $f(\boldsymbol{v})=\phi(\boldsymbol{A}\boldsymbol{v})$, for some function $\phi:\mathbb{R}^{s}\rightarrow\mathbb{R}_{+}$ and matrix $\boldsymbol{A}\in\mathbb{R}^{s\times n}$.

To identify the analysis/synthesis dichotomy in terms of how the patches are related to the underlying image, not in terms of how the patches are modelled, we refer to the formulation reviewed in this section as patch-analysis and to its synthesis counterpart (to be introduced in the next section) as patch-synthesis.

4 Patch-Synthesis Formulation

We now present the patch-synthesis formulation, which can be summarized as follows: the clean image is generated by additively combining a collection of patches; the patches themselves follow some probabilistic model, but are a priori mutually independent.

Consider a collection of patches $\{\boldsymbol{z}_{m}\in\mathbb{R}^{n},\;m=1,...,M\}$ and let image $\boldsymbol{x}$ be synthesized from these patches by combining them additively according to

[TABLE]

where matrices $\boldsymbol{Q}_{m}\in\mathbb{R}^{N\times n}$ are such that they average the values in the several patches that contribute to a given pixel of $\boldsymbol{x}$.

A simple 1D example will help clarify the structure of the $\boldsymbol{Q}_{m}$ matrices. Consider that $\boldsymbol{x}=[x_{1},x_{2},x_{3},x_{4}]^{T}$ is produced by combining all the consecutive 2-element patches (with periodic boundary conditions), which correspond to the subsets of components $\{1,2\}$, $\{2,3\}$, $\{3,4\},$ and $\{4,1\}$. That is, $M=4$ and

[TABLE]

Stacking all the patches in vector $\boldsymbol{z}\in\mathbb{R}^{Mn}$ and considering a matrix $\boldsymbol{Q}=[\boldsymbol{Q}_{1}\cdots\boldsymbol{Q}_{M}]\in\mathbb{R}^{N\times(Mn)}$, the synthesis expression in (15) can be written compactly as

[TABLE]

With the patches modelled as independent and identically distributed samples of some patch-wise pdf $g$, their joint log-prior is

[TABLE]

Notice that obtaining samples from (19) is as simple as obtaining samples from $g$ itself. Consequently, generating image samples under this synthesis model simply corresponds to generating samples from (19) and then multiplying them by $\boldsymbol{Q}$. This is in contrast with the patch-analysis prior (6), where, even if $f$ is a valid pdf, it is not trivial to obtain samples from $p_{X}$.

The resulting MAP denoising criterion can now be written as

[TABLE]

As the patch-analysis model, the patch-synthesis formulation that we have just presented is compatible with both analysis and synthesis priors for the patches, and of course with any valid pdf for vectors in $\mathbb{R}^{n}$. An analysis formulation simply amounts to choosing a patch prior of the form $g=\psi\circ\boldsymbol{B}$, i.e., $g(\boldsymbol{z}_{m})=\psi(\boldsymbol{B}\boldsymbol{z}_{m})$, for some function $\psi:\mathbb{R}^{s}\rightarrow\mathbb{R}_{+}$ and matrix $\boldsymbol{B}\in\mathbb{R}^{s\times n}$.

In a synthesis formulation, each patch $\boldsymbol{z}_{m}$ is synthesized using some dictionary $\boldsymbol{D}$ as \boldsymbol{z}_{m}=\boldsymbol{D}\mbox{\boldmath\gamma}_{m}, and the \mbox{\boldmath\gamma}_{m} follow some prior $p_{\Gamma}$; stacking all the patch coefficients in vector $\gamma$, we can write \boldsymbol{z}=\mbox{\boldmath\Delta}\mbox{\boldmath\gamma} (where $\Delta$ is a block-diagonal matrix with $M$ replicas of $\boldsymbol{D}$), thus the MAP denoising criterion becomes

[TABLE]

where

[TABLE]

5 Relationship Between the Analysis and Synthesis Formulations

Leaving aside for now the choice of the patch priors, let us focus on the relationship between formulations (13) and (20). The key observation underlying the relationship between these two formulations is that (assuming the patch structure in both formulations is the same)

[TABLE]

where $\boldsymbol{I}$ denotes the identity matrix, but in general

[TABLE]

In other words, $\boldsymbol{Q}$ is a left pseudo-inverse of $\boldsymbol{P}$. To prove (23), simply notice that if a collection of patches is extracted from some image and then these patches are used to synthesize an image by averaging the overlapping pixels, an identical image is obtained. Of course, the converse is not true, in general: if an image is synthesized by averaging the overlapping pixels of a collection of patches, and then patches are extracted from the synthesized image, there is no guarantee that these patches are equal to the original ones, thus proving (24). In the trivial and uninteresting cases where the patches are singletons, or non-overlapping, we would have $\boldsymbol{P}\,\boldsymbol{Q}=\boldsymbol{I}$.

As shown in [7], given the analysis formulation (13), an equivalent synthesis formulation is

[TABLE]

where the constraint $\boldsymbol{z}\in\mathop{\mbox{range}}(\boldsymbol{P})$ enforces $\boldsymbol{z}$ to be in the subspace spanned by the columns of $\boldsymbol{P}$. Notice that this constraint corresponds to having a collection of patches extracted from some image, i.e., it forces the patches to agree on the value of each shared pixel. Since (20) does not enforce this constraint, it is not, in general, equivalent to the patch-analysis formulation (13).

6 ADMM for Patch-Synthesis Denoising

In this section, we derive an instance of ADMM to deal with (20). The first step is to rewrite it as a constrained problem,

[TABLE]

where

[TABLE]

is the negative log-prior, or regularizer. ADMM for this problem takes the form

[TABLE]

The update equation for $\boldsymbol{u}^{(t+1)}$ is a denoising step. Due to the separability of the squared $\ell_{2}$ norm and of $\xi$ (see (27)), this can be separately solved with respect to each patch: for $m=1,...,M$,

[TABLE]

which is simply the proximity operator of $(1/\rho)\xi_{m}$ [1], computed at $\boldsymbol{d}_{m}^{(t)}-\boldsymbol{z}_{m}^{(t+1)}$. In a Bayesian viewpoint, (28) corresponds to obtaining the MAP estimate of $\boldsymbol{u}_{m}$ from observations $\boldsymbol{d}_{m}^{(t)}-\boldsymbol{z}_{m}^{(t+1)}$, assuming additive white Gaussian noise of variance $1/\rho$ and a negative log-prior $\xi_{m}$. If the $\xi_{m}$ are convex, this MAP estimate is unique, due to the strict convexity of the quadratic term in (28).

Computing $\boldsymbol{z}^{(t+1)}$ corresponds to solving an unconstrained quadratic problem, the solution being

[TABLE]

The bottleneck in this update equation seems to be the matrix inversion, since $\boldsymbol{Q}^{T}\boldsymbol{Q}$ is a huge $(Mn)\times(Mn)$ matrix. However, this inversion can be solved very efficiently by resorting to the Sherman-Morrison-Woodbury matrix inversion formula. In fact,

[TABLE]

where matrix $(\sigma^{2}\rho\,\boldsymbol{I}+\boldsymbol{Q}\,\boldsymbol{Q}^{T})$ is diagonal, thus its inversion is trivial. To prove that $\boldsymbol{Q}\boldsymbol{Q}^{T}$ is diagonal, recall that $\boldsymbol{Q}=\begin{bmatrix}\boldsymbol{Q}_{1}\;\boldsymbol{Q}_{2}\;\cdots\boldsymbol{Q}_{M}\end{bmatrix}$, thus

[TABLE]

The element $(i,j)$ of matrix $\boldsymbol{Q}_{m}\boldsymbol{Q}_{m}^{T}\in\mathbb{R}^{N\times N}$ is the inner product between the $i$-th and the $j$-th rows of $\boldsymbol{Q}_{m}$. Since each pixel in the $m$-th patch contributes to one and only one pixel in the synthesized image, the rows of $\boldsymbol{Q}_{m}$ have disjoint support, thus $(i\neq j)\Rightarrow(\boldsymbol{Q}_{m}\boldsymbol{Q}_{m}^{T})_{i,j}=0$, that is, $\boldsymbol{Q}_{m}\boldsymbol{Q}_{m}^{T}$ has no non-zero elements outside of its main diagonal, thus is a diagonal matrix. Finally, since $\boldsymbol{Q}\boldsymbol{Q}^{T}$ is a sum of diagonal matrices, it is a diagonal matrix.

It is also easy to obtain explicitly the elements of the diagonal of $\boldsymbol{Q}\boldsymbol{Q}^{T}$. The diagonal elements $(\boldsymbol{Q}_{m}\boldsymbol{Q}_{m}^{T})_{i,i}$ are given by

[TABLE]

If the $m$-th patch contributes to pixel $i$, the sum in (32) contains exactly one non-zero term, thus it is equal to the square of the weight with which the $m$-th patch contributes to the synthesis of pixel $i$. If the $m$-th patch does not contributes to pixel $i$, then $(\boldsymbol{Q}_{m}\boldsymbol{Q}_{m}^{T})_{i,i}=0$. Finally, since the weight with which each patch element contributes to each pixel equals the inverse of the number of patches that contributes to that pixel, $(\boldsymbol{Q}\boldsymbol{Q}^{T})_{i,i}$ is equal to the inverse of the number of patches that contribute to pixel $i$. Referring to the example in (16)–(17), we have simply $\boldsymbol{Q}\boldsymbol{Q}^{T}=(1/2)\boldsymbol{I}$, because each pixel is synthesized from two patches.

Finally, letting \boldsymbol{q}=\mbox{diag}\bigl{(}(\sigma^{2}\rho\boldsymbol{I}+\boldsymbol{Q}\boldsymbol{Q}^{T})^{-1}\bigr{)}, and denoting $\boldsymbol{s}^{(t)}=\boldsymbol{Q}^{T}\boldsymbol{y}+\sigma^{2}\rho(\boldsymbol{u}^{(t)}+\boldsymbol{d}^{(t)})$, we can write the update equation (29) as

[TABLE]

where $\odot$ denotes element-wise product between two vectors. The leading cost of this update is that of the matrix-vector products involving $\boldsymbol{Q}$ and $\boldsymbol{Q}^{T}$, which is $O(MNn)$; all the other operations in (29) have lower computational cost.

7 Conclusions and Future Work

In this paper, we have revisited patch-based image priors under the light of the synthesis vs analysis dichotomy. After showing that the classical patch-based image models (namely the EPLL) corresponds to an analysis formulation, we have proposed a patch-synthesis formulation, and analyzed its relationship with the analysis formulation, showing that they are, in general, not equivalent. Finally, we have shown how to address image denoising under the proposed formulation, via an ADMM algorithm.

We stress again that the purpose of this paper is not to introduce a new particular image prior, but a general patch-based synthesis formulation/framework, which (to the best of our knowledge) was missing from the literature on patch-based image processing, and which can be instantiated with many different patch models/priors. For this reason, we have abstained from presenting experimental results; these would critically depend on the choice and estimation of a particular patch model, which is not the focus of this paper.

Ongoing work includes the development of efficient algorithms for learning patch models under the assumption that they will be used in a synthesis formulation.