Scene-adapted plug-and-play algorithm with convergence guarantees

TL;DR

This paper introduces a scene-adapted plug-and-play algorithm using a Gaussian mixture model denoiser within an ADMM framework, providing convergence guarantees and applying it to hyperspectral sharpening.

Contribution

It develops a convergent plug-and-play method with scene-specific priors, advancing the application of denoising models in image processing tasks.

Findings

Proves convergence of the proposed algorithm.

Demonstrates effectiveness in hyperspectral sharpening.

Utilizes scene-adapted priors for improved performance.

Abstract

Recent frameworks, such as the so-called plug-and-play, allow us to leverage the developments in image denoising to tackle other, and more involved, problems in image processing. As the name suggests, state-of-the-art denoisers are plugged into an iterative algorithm that alternates between a denoising step and the inversion of the observation operator. While these tools offer flexibility, the convergence of the resulting algorithm may be difficult to analyse. In this paper, we plug a state-of-the-art denoiser, based on a Gaussian mixture model, in the iterations of an alternating direction method of multipliers and prove the algorithm is guaranteed to converge. Moreover, we build upon the concept of scene-adapted priors where we learn a model targeted to a specific scene being imaged, and apply the proposed method to address the hyperspectral sharpening problem.

Tables

Table 1. Table 1 : HS and MS fusion on (cropped) ROSIS Pavia University and Moffett Field datasets.

		Exp. 1 (PAN)			Exp. 2 (PAN)			Exp. 3 (R,G,B,N-IR)			Exp. 4 (R,G,B,N-IR)
Dataset	Metric	ERGAS	SAM	SRE	ERGAS	SAM	SRE	ERGAS	SAM	SRE	ERGAS	SAM	SRE
Rosis	Dictionary [14]	1.99	3.28	22.64	2.05	3.16	22.32	0.47	0.85	34.60	0.85	1.47	29.66
	GMM [15]	1.75	2.89	23.67	1.92	2.92	22.85	0.48	0.87	34.32	0.91	1.65	29.05
	proposed GMM	1.65	2.75	24.17	1.81	2.76	23.31	0.49	0.87	34.59	0.80	1.42	30.14
Moffett	Dictionary [14]	2.67	4.18	20.28	2.74	4.20	20.05	1.85	2.72	23.58	2.12	3.21	22.25
	GMM [15]	2.66	4.24	20.26	2.78	4.27	19.87	1.81	2.68	23.81	1.98	2.93	22.91
	proposed GMM	2.54	4.06	20.66	2.65	4.10	20.28	1.73	2.58	24.18	1.97	2.90	22.94

Equations

$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+g(\mathbf{H}\mathbf{x}),$$

$${\mathbf{x}}^{(k+1)}$$

$${\mathbf{v}}^{(k+1)}$$

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

$$\underset{\mathbf{x}}{\min }\sum _{j=1}^J{g}_j({\mathbf{H}}_j\mathbf{x})$$

$$f(\mathbf{x})=0,\mathbf{H}=\left[\begin{array}{c}{\mathbf{H}}_1\\ ⋮\\ {\mathbf{H}}_J\end{array}\right],\mathbf{v}=\left[\begin{array}{c}{\mathbf{v}}_1\\ ⋮\\ {\mathbf{v}}_J\end{array}\right],g(\mathbf{v})=\sum _{j=1}^J{g}_j({\mathbf{v}}_j).$$

$${\mathbf{x}}^{(k+1)}$$

$${\mathbf{v}}_1^{(k+1)}$$

$$⋮$$

$${\mathbf{v}}_J^{(k+1)}$$

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

$$p({\bf x}_{i})=\sum_{j=1}^{K}\alpha_{j}\mathcal{N}({\bf x}_{i};\mbox{\boldmath$\mu$}_{j},{\bf C}_{j}),$$

$${\hat{\mathbf{x}}}_i=\sum _{j=1}^K{\beta }_j({\mathbf{y}}_i){\mathbf{v}}_j({\mathbf{y}}_i),$$

$${\bf v}_{j}({\bf y}_{i})={\bf C}_{j}\Bigl{(}{\bf C}_{j}+\sigma^{2}{\bf I}\Bigr{)}^{-1}{\bf y}_{i},$$

$${\beta }_j({\mathbf{y}}_i)=\frac{{\alpha }_j\mathcal{N}({\mathbf{y}}_i;0,{\mathbf{C}}_j+{\sigma }^2\mathbf{I})}{{\sum }_{k=1}^K{\alpha }_j\mathcal{N}({\mathbf{y}}_i;0,{\mathbf{C}}_k+{\sigma }^2\mathbf{I})}.$$

$$\hat{\mathbf{x}}\in \underset{\mathbf{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}\sum _{i=1}^N\Vert {\hat{\mathbf{x}}}_i-{\mathbf{P}}_i\mathbf{x}{\Vert }_2^2,$$

$$\hat{{\bf x}}=\Bigl{(}\sum_{i=1}^{N}{\bf P}_{i}^{T}{\bf P}_{i}\Bigr{)}^{-1}\Bigl{(}\sum_{i=1}^{N}{\bf P}_{i}^{T}\hat{\bf x}_{i}\Bigr{)}=\frac{1}{n_{p}}\,\sum_{i=1}^{N}{\bf P}_{i}^{T}\hat{\bf x}_{i},$$

$${\hat{\mathbf{x}}}_i$$

$$\hat{\mathbf{x}}=\underset{\mathbf{W}}{\underbrace{\frac{1}{n_p}\sum _{i=1}^N{\mathbf{P}}_i^T{\mathbf{F}}_i{\mathbf{P}}_i}}\mathbf{y}=\mathbf{W}\mathbf{y}.$$

$$\text{mbox}{prox}_g(\mathbf{y})=\arg \underset{\mathbf{x}}{\min }\frac{1}{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2^2+\frac{1}{2}{\mathbf{x}}^T({\mathbf{W}}^{-1}-\mathbf{I})\mathbf{x}.$$

$$\frac{1}{2}\,{\bf x}^{T}\bigl{(}{\bf I}+{\bf W}^{-1}-{\bf I}\bigr{)}{\bf x}-{\bf x}^{T}{\bf y}+\frac{1}{2}\,\|{\bf y}\|_{2}^{2},$$

$$\mathbf{H}=\left[\begin{array}{c}{\mathbf{H}}_1\\ ⋮\\ {\mathbf{H}}_{J-1}\\ \mathbf{I}\end{array}\right],$$

$$\Vert f(\mathbf{x})-f(\mathbf{y}){\Vert }^2\le \Vert \mathbf{x}-\mathbf{y}{\Vert }^2.$$

$${\mathbf{Y}}_h=\mathbf{Z}\mathbf{B}\mathbf{M}+{\mathbf{N}}_h,{\mathbf{Y}}_m=\mathbf{R}\mathbf{Z}+{\mathbf{N}}_m,$$

$${\mathbf{Y}}_h=\mathbf{E}\mathbf{X}\mathbf{B}\mathbf{M}+{\mathbf{N}}_h,{\mathbf{Y}}_m=\mathbf{R}\mathbf{E}\mathbf{X}+{\mathbf{N}}_m.$$

$$\hat{\mathbf{X}}\in \underset{\mathbf{X}}{\mathrm{argmin}}\frac{\Vert \mathbf{E}\mathbf{X}\mathbf{B}\mathbf{M}-{\mathbf{Y}}_h{\Vert }_F^2}{2}+\frac{\lambda \Vert \mathbf{R}\mathbf{E}\mathbf{X}-{\mathbf{Y}}_m{\Vert }_F^2}{2}+\tau \varphi (\mathbf{X}),$$

$$\Vert \mathbf{E}\mathbf{X}\mathbf{B}\mathbf{M}-{\mathbf{Y}}_h{\Vert }_F^2$$

$$\Vert \mathbf{R}\mathbf{E}\mathbf{X}-{\mathbf{Y}}_m{\Vert }_F^2$$

$$g_1({\mathbf{v}}_1)$$

$$g_2({\mathbf{v}}_2)$$

$$g_3({\mathbf{v}}_3)$$

$${\mathbf{H}}_1=({\mathbf{B}}^T\otimes \mathbf{I}),{\mathbf{H}}_2=\mathbf{I},{\mathbf{H}}_3=\mathbf{I}.$$

$${\mathbf{X}}^{k+1}=\arg \underset{\mathbf{X}}{\min }\Vert \mathbf{X}\mathbf{B}-{\mathbf{V}}_1^k-{\mathbf{D}}_1^k{\Vert }_F^2+\Vert \mathbf{X}-{\mathbf{V}}_2^k-{\mathbf{D}}_2^k{\Vert }_F^2$$

$$+\Vert \mathbf{X}-{\mathbf{V}}_3^k-{\mathbf{D}}_3^k{\Vert }_F^2,$$

$${\mathbf{V}}_1^{k+1}=\arg \underset{{\mathbf{V}}_1}{\min }\Vert \mathbf{E}{\mathbf{V}}_1\mathbf{M}-{\mathbf{Y}}_h{\Vert }_F^2+\rho \Vert {\mathbf{X}}^{k+1}\mathbf{B}-{\mathbf{V}}_1-{\mathbf{D}}_1^k{\Vert }_F^2,$$

$${\mathbf{V}}_2^{k+1}=\arg \underset{{\mathbf{V}}_2}{\min }\lambda \Vert \mathbf{R}\mathbf{E}{\mathbf{V}}_2-{\mathbf{Y}}_m{\Vert }_F^2+\rho \Vert {\mathbf{X}}^{k+1}-{\mathbf{V}}_2-{\mathbf{D}}_2^k{\Vert }_F^2,$$

Full text

Index Terms—  Plug-and-play, Gaussian mixture denoiser, scene-adapted prior, data fusion, hyperspectral sharpening.

1 Introduction

Image denoising is not only one of the core problems in image processing, but also a building block for many other tasks, and has been a very active research topic over the last decades. Most of the current state-of-the-art methods are patch-based, i.e., they rely on a divide and conquer principle, where instead of dealing with the noisy image as a whole, small (overlapping) patches are extracted, denoised independently, and put back in their locations [1, 2, 3]. Applying these patch-based methods in more general restoration/reconstruction problems is not trivial and has been a topic of research in the past years [4, 5, 6].

Recently, Venkatakrishnan et al. proposed a flexible framework, called plug-and-play (PnP), where a state-of-the-art denoiser is seen as a black box and plugged into the iterations of an alternating direction method of multipliers (ADMM) [7]. This approach allows using an arbitrary denoiser as a regularizer in an imaging inverse problem, such as deblurring or super-resolution, tackled via an ADMM algorithm, but departing from its standard use, where the denoiser is the proximity operator of some convex regularizer [8]. However, plugging an arbitrary denoiser (possibly without a closed-form expression) into ADMM begs obvious questions [9, 7]: is the resulting algorithm guaranteed to converge? If so, does it converge to a (global or local) optimum of some objective function? Can we identify this function? Here, we give positive answers to these questions, when the plugged denoiser is a modified Gaussian mixture model (GMM) based denoiser [2, 3, 10].

As proposed in earlier work [11], the GMM can be adapted to specific classes of images; the rationale is that denoisers based on such class-adapted priors are able to better capture the characteristics of the class than a general-purpose denoiser. In this paper, we take this adaptation one step further, by considering scene-adapted priors. As the name suggests, the model is no longer learned from a set of images from the same class as the input image, but from one or more images over the same scene being restored/reconstructed. Hyperspectral sharpening [12, 13, 14, 15] and image deblurring with noisy/blurred image pairs [16] are two examples of applications that may leverage from such priors.

In summary, the contributions are threefold: (i) we prove that the ADMM algorithm with a plugged-in GMM-based denoiser is guaranteed to converge to the global minimum of a cost function, provided that we introduce a small modification to the denoiser; (ii) this proposed modification not only guarantees convergence of the algorithm, but it is also motivated by the scene-adapted perspective; (iii) we show that the proposed algorithm yields results that are often better than other state-of-the-art hyperspectral sharpening methods.

The paper is organized as follows. Section 2 gives a brief description of ADMM, which is the basis of the proposed algorithm. Section 3 explains the rationale behind the PnP framework, and our reasons for choosing the GMM-based denoiser. Section 4 contains what is arguably the main contribution of this paper: the convergence of the resulting algorithm. The application to HS sharpening is described in Section 5, and experimental results are reported in Section 6. Finally, Section 7 concludes the paper.

2 Alternating Direction Method of Multipliers

Although dating back to the 1970’s [17], ADMM has seen a surge of interest in the last decade, as a flexible and efficient optimization tool, widely used in imaging problems, machine learning, and other areas [18, 19]. One of the canonical problems for ADMM has the form

[TABLE]

where $f$ and $g$ are closed, proper, convex functions, and matrix ${\bf H}$ is of appropriate dimensions [18]. Each iteration of ADMM for (1) is as follows (with the superscript $(\cdot)^{(k)}$ denoting the iteration counter):

[TABLE]

where ${\bf u}^{(k)}$ are the (scaled) Lagrange multipliers at iteration $k$, and $\rho$ is the penalty parameter.

ADMM can be directly applied to problems involving the sum of $J$ closed, proper, convex functions, composed with linear operators,

[TABLE]

by casting it into the form (1) as follows:

[TABLE]

The resulting instance of ADMM (called SALSA–split augmented Lagrangian shrinkage algorithm [20]) is

[TABLE]

where (as for ${\bf v}^{(k)}$), {\bf u}^{(k)}=\bigl{[}\bigl{(}{\bf u}_{1}^{(k)}\bigr{)}^{T},\,\dots,\,\bigl{(}{\bf u}_{J}^{(k)}\bigr{)}^{T}\bigr{]}^{T}.

3 Plug-and-Play Priors

The PnP framework [7] emerges from noticing that sub-problems (7)-(8) can be interpreted as the MAP solutions to pure denoising problems. In such situations, the resulting sub-problem has an identity observation operator, function $g_{i}$ acting as the regularizer (or negative log-prior, in Bayesian terms), $({\bf x}^{(k+1)}\!-{\bf u}^{(k)})$ as the noisy data, and noise variance equal to $1/\rho$. The goal is to capitalize on recent developments in image denoising, by using a state-of-the-art denoiser plugged into the iterations of ADMM, instead of trying to come up with a convex regularizer for which it is possible to compute the corresponding proximity operator.

3.1 GMM-Based Denoising

Amongst the several state-of-the-art denoising methods, we adopt a GMM patch-based denoiser. The reasons for this choice are threefold: (i) it has been shown that GMM are good priors for clean image patches [3, 10]; (ii) a GMM can be learned either from an external dataset of clean images [3, 10], or directly from noisy patches [2]; (iii) the ease of learning a GMM from an external dataset opens the door to class-specific priors, which, naturally, lead to superior results since they capture better the characteristics of the image class in hand, than a general-purpose image model [11, 21].

Leaving aside for now the question of how to learn the parameters of the model, we consider a GMM prior111\mathcal{N}({\bf x};\mbox{\boldmath\mu},{\bf C}) denotes a Gaussian probability density function of mean $\mu$ and covariance $\bf C$, computed at ${\bf x}$.

[TABLE]

for any patch ${\bf x}_{i}$ of the image to be denoised. A current practice in patch-based image denoising is to treat the average of each patch separately: the average of the noisy patch is subtracted from it, the resulting zero-mean patch is denoised, and the average is added back to the estimated patch. The implication of this procedure is that, without loss of generality, we can assume that the components of the GMM have zero mean (\mbox{\boldmath\mu}_{j}=0, for $j=1,...,K$).

Given this GMM prior and a noisy version ${\bf y}_{i}$ of each patch, where ${\bf y}_{i}={\bf x}_{i}+{\bf n}_{i}$, with ${\bf n}_{i}\sim\mathcal{N}(0,\sigma^{2})$, the minimum mean squared error (MMSE) estimate of ${\bf x}_{i}$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

Notice that $\beta_{j}({\bf y}_{i})$ is the posterior probability that the i-th patch was generated by the j-th component of the GMM, and ${\bf v}_{j}({\bf y}_{i})$ is the MMSE estimate of the i-th patch if it was known that it had been generated by the $j$-th component.

As is standard in patch-based denoising, after computing the MMSE estimates of the patches, they are returned to their location and combined by straight averaging. This corresponds to solving the following optimization problem

[TABLE]

where ${\bf P}_{i}\in\{0,1\}^{n_{p}\times n}$ is a binary matrix that extracts the $i$-th patch from the image (thus ${\bf P}_{i}^{T}$ puts the patch back into its place), $N$ is the number of patches, $n_{p}$ is the number of pixels in each patch, and $n$ the total number of image pixels. The solution to (14) is

[TABLE]

assuming the patches are extracted with unit stride and periodic boundary conditions, thus every pixel belongs to $n_{p}$ patches and $\sum_{i=1}^{N}{\bf P}_{i}^{T}{\bf P}_{i}=n_{p}{\bf I}$.

3.2 From Class-Adapted to Scene-Adapted Priors

Class-adapted priors have been used successfully in image deblurring, if the image to be deblurred is known to belong to a certain class, e.g., image of a face or fingerprint, and a prior adapted to that class is used [11]. Here, we propose taking this adaptation to an extreme, whenever we have two types of data describing the same scene, as is the case with data fusion problems. Instead of leveraging a dataset of images from the same class as the input image, we should leverage one type of data to learn a model that is targeted to the specific scene being imaged, and use it as a prior in the reconstruction/restoration process. The underlying assumption is that the two types of data share the same spatial statistical properties; this is a reasonable assumption, since the scene being imaged is exactly the same. In Section 6, we implement this idea to tackle a data fusion problem known as hyperspectral (HS) sharpening (see [12] for a detailed description of the problem). We resort to a GMM prior, as described in Subsection 3.1, learned from (patches of the bands of) the observed multispectral (MS) or panchromatic (PAN) image using the classical expectation-maximization (EM) algorithm, and then leverage this prior to sharpen the low resolution HS bands.

In light of the scene-adaptation scheme just described, we further claim that, it not only makes sense to use a GMM prior learned from the observed MS or PAN image, but it also makes sense to keep the weights $\beta_{m}$ of each patch in the target image, used to compute (11), equal to its value for the corresponding observed patch, obtained during training. In other words, we use the weights as if we were simply denoising the high resolution MS image, with a GMM trained from its noisy patches. In fact, another instance of this idea has been previously used, based on sparse representations on learned dictionaries [14]. We stress that this approach is only reasonable in the case of image denoising or image reconstruction using scene-adapted priors, where the training and target images are of the same scene, and thus have similar spatial structures. As will be shown in the next section, this modification will play an important role in proving convergence of the proposed PnP-ADMM algorithm.

4 Convergence of the PnP-ADMM Algorithm

We begin by presenting a simplified version of a classical theorem on the convergence of a generalized version of ADMM, proved in the seminal paper by Eckstein and Bertsekas [18]. The full version of the theorem allows for inexact solution of the optimization problems in the iterations of ADMM, while this version, based on which our proof of convergence of PnP-ADMM will be supported, assumes exact solutions.

Theorem 1** (Eckstein and Bertsekas [18]).**

Consider a problem of the form (1), where ${\bf H}$ has full column rank, and $f:\mathbb{R}^{n}\rightarrow\bar{\mathbb{R}}$ and $g:\mathbb{R}^{m}\rightarrow\bar{\mathbb{R}}$ are closed, proper, and convex functions, and let ${\bf v}_{0},{\bf u}_{0}\in\mathbb{R}^{m}$, and $\rho>0$ be given. If the sequences $({\bf x}^{(k)},k=0,1,\dots)$, $({\bf v}^{(k)},k=0,1,\dots)$, and $({\bf u}^{(k)},k=0,1,\dots)$ are generated according to (2), (3), (4), then $({\bf x}^{(k)},k=0,1,\dots)$ converges to a solution of (1), ${\bf x}^{(k)}\rightarrow{\bf x}^{*}$, if one exists. Furthermore, if a solution does not exist, then at least one of the sequences $({\bf v}^{(k)},k=0,1,\dots)$ or $({\bf u}^{(k)},k=0,1,\dots)$ diverges.

Before proceeding, some properties of the denoising function are established.

4.1 Analysis of GMM-Based Denoiser

Consider the role of the denoiser: for a noisy input argument ${\bf y}\in\mathbb{R}^{n}$, it produces an estimate $\hat{{\bf x}}\in\mathbb{R}^{n}$, given the parameters of the GMM, $\theta=\{{\bf C}_{1},\dots,{\bf C}_{K},\beta_{1}^{1},\dots,\beta_{K}^{n}\}$, and a noise level, $\sigma$. As the GMM-denoiser is patch-based, the algorithm starts by extracting (overlapping, with unit stride) patches, followed by computing their estimates using (11). Recall that, with the scene-adaptation scheme, parameters $\beta$ no longer depend on the noisy patches, and thus we may write compactly, for each patch,

[TABLE]

where ${\bf P}_{i}\in\{0,1\}^{n_{p}\times n}$ is the binary matrix defined in Subsection 3.1, thus ${\bf y}_{i}={\bf P}_{i}\,{\bf y}$. Moreover, we can relate the input of the denoiser to its image-wise output using (15), under the assumption of periodic boundary conditions. Once again, using compact notation

[TABLE]

which shows that the GMM-based MMSE denoiser, with fixed weights, is a linear function of the noisy image. Matrix ${\bf W}$, of course, depends on the parameters of the method ($\sigma$, $n_{p}$) as well as on the parameters of the GMM, but it is fixed throughout the iterations of the ADMM loop. Recall that the weights, $\beta$, and covariance matrices, ${\bf C}$, are obtained during training, and $\sigma$ is assumed to be known, thus is fixed. The next lemma states the key properties of matrix ${\bf W}$. The proof can be found in the Appendix.

Lemma 1**.**

Consider that all covariance matrices $\{{\bf C}_{1},...,{\bf C}_{K}\}$ in the GMM are positive definite. Then, ${\bf W}$ is symmetric, positive definite, and has spectral norm strictly less than 1.

Because the denoiser corresponds to multiplying the argument by a symmetric positive definite matrix, it is the proximity operator of a quadratic function, as stated in the following lemma.

Lemma 2**.**

The function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, defined by $f({\bf y})={\bf W}{\bf y}$, where ${\bf W}$ is symmetric, positive definite, and has spectral norm strictly less than 1, corresponds to the proximity operator of the quadratic funtion $g({\bf y})=\frac{1}{2}{\bf y}^{T}({\bf W}^{-1}-{\bf I}){\bf y}$.

Proof: Using the definition of proximity operator [22],

[TABLE]

The function being minimized w.r.t. $\bf x$ can be written as

[TABLE]

which is quadratic, strictly convex, with a unique minimizer ${\bf W}{\bf y}$. $\blacksquare$

4.2 Convergence Proof

Lemmas 1 and 2 allow proving convergence of the proposed algorithm, which is formalized in the following theorem.

Theorem 2**.**

Consider the application of SALSA to (5), which results in (6) to (9). Consider also that ${\bf v}_{1},\dots,{\bf v}_{J-1}$ are proximity operators, ${\bf v}_{J}$ is obtained with the GMM-denoiser described in 3.1, and that all the covariance matrices $\{{\bf C}_{1},...,{\bf C}_{K}\}$ in the GMM are positive definite. Then, the proposed algorithm converges.

Proof: The proof consists in verifying the conditions of Theorem 1.

The first condition for convergence is that ${\bf H}$ has full column rank. The proposed algorithm is an instance of SALSA (as defined in (6)–(9)), which in turn is an instance of ADMM (as defined in (2)–(4)), with

[TABLE]

which has full column rank, regardless of ${\bf H}_{1},\dots,{\bf H}_{J-1}$.

The second condition regards functions $g_{1}$ to $g_{J}$. In fact, $g_{1}$ to $g_{J-1}$ are closed, proper, and convex, by virtue of ${\bf v}_{1}$ to ${\bf v}_{J-1}$ being proximity operators [22]. Furthermore, as shown by Lemma 2, the denoiser also corresponds to the proximity operator of a quadratic (thus also closed and proper) function $g_{J}$, which is strongly convex. Finally, the function being minimized is the sum of $J$ convex functions, one of which is strongly convex, thus a minimizer is guaranteed to exist. $\blacksquare$

We conclude by observing that this proof of convergence depends critically on the fact that the denoiser uses fixed weights in $\beta$, rather than being a pure MMSE denoiser, as described in Subsection 3.1. In fact, the simple univariate example in Fig. 1a shows that an MMSE denoiser may not even be non-expansive222Recall that a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is non-expansive if, for any ${\bf x},{\bf y}\in\mathbb{R}^{n}$,

$\|f({\bf x})-f({\bf y})\|^{2}\leq\|{\bf x}-{\bf y}\|^{2}.$

(18)

An important connection between non-expansiveness and proximity operators in given in the following theorem:

Theorem 3** (Moreau [23]).**

A function is a proximal mapping if and only if it is non-expansive and the sub-gradient of some convex function.

. This figure plots the MMSE estimate $\hat{x}$ in (11), as a function of the noisy observed $y$, under a GMM prior with two zero-mean components. Clearly, the presence of regions of the function with derivative larger than 1 shows that it is not non-expansive.

5 Application: Hyperspectral Sharpening

In this section, we tackle a data fusion problem known as HS sharpening, under the PnP framework with a GMM-denoiser and a scene-adapted prior, as described above. Using compact matrix notation (as in [13]), the HS sharpening inverse problem can be modeled as

[TABLE]

where ${\bf Z}\in\mathbb{R}^{L_{h}\times n_{m}}$ is the target image to be estimated, ${\bf Y}_{h}\in\mathbb{R}^{L_{h}\times n_{h}}$ is the observed HS data, ${\bf B}\in\mathbb{R}^{n_{m}\times n_{m}}$ is a spatial convolution operator, ${\bf M}\in\mathbb{R}^{n_{m}\times n_{h}}$ is a sub-sampling operator, ${\bf Y}_{m}\in\mathbb{R}^{L_{m}\times n_{m}}$ is the observed MS data, ${\bf R}\in\mathbb{R}^{L_{m}\times L_{h}}$ models the spectral responses of the MS sensor, and ${\bf N}_{h}$ and ${\bf N}_{m}$ are Gaussian noises of known variances. Moreover, we assume the sensor has a fixed point spread function, and matrices ${\bf B}$ and ${\bf R}$ are known; see [13] for a blind approach.

HS data typically has hundreds of bands but the spectral vectors, i.e. the columns of ${\bf Z}$, live in a low dimensional subspace which may be learned from the observed HS data [24]. In particular, for independent and identically distributed (iid) noise, the range of ${\bf Z}$ may be learned from the eigenvectors of the correlation matrix ${\bf Y}_{h}{\bf Y}_{h}^{T}/n_{h}$. Then, instead of estimating ${\bf Z}$ directly, we estimate latent images ${\bf X}\in\mathbb{R}^{L_{s}\times n_{m}}$, and recover ${\bf Z}={\bf E}{\bf X}$, where the columns of ${\bf E}\in\mathbb{R}^{L_{h}\times L_{s}}$ are the eigenvectors corresponding to the $L_{s}$ largest eigenvalues of ${\bf Y}_{h}{\bf Y}_{h}^{T}/n_{h}$. The observation model becomes

[TABLE]

A classical approach to this inverse problem is to compute the maximum a posteriori (MAP) estimate, i.e.,

[TABLE]

where $\phi$ is the negative log-prior (or regularizer), while $\lambda$ and $\tau$ control the relative weight of each term. We follow Simões et al. [13] and use SALSA to tackle (21), where once again we assume the noise is Gaussian iid; more general forms of noise correlation can be included [14]. Notice that, while in (21) the optimization variable is a matrix ${\bf X}$, in (5) the optimization variable is a vector ${\bf x}$. These is merely a difference in notation, as we can represent matrix ${\bf X}$ by its vectorized version ${\bf x}=\mbox{vec}({\bf X})$, which corresponds to stacking the columns of ${\bf X}$ into a vector. Using well-known equalities relating vectorization and Kronecker products, we can write

[TABLE]

and map (21) into (5) by letting $J=3$, and

[TABLE]

where ${\bf v}_{i}=\mbox{vec}({\bf V}_{i})$, for $i=1,2,3$, while

[TABLE]

Although we can arbitrarily switch between matrix and vector representations, in what follows it is more convenient to use matrix notation. Each iteration of the resulting ADMM has the form

[TABLE]

where $\rho$ is the so-called penalty parameter of ADMM, and ${\bf D}_{i}$ are the scaled Lagrange dual variables [19]. The first three problems are quadratic, thus with closed-form solutions involving matrix inversions (for details, see [13]):

[TABLE]

where $\odot$ denotes entry-wise (Hadamard) product. The matrix inversion in (27) can be implemented with cost $O(n_{m}\log n_{m})$ using the FFT, under the assumption of periodic boundary conditions. The matrix inversions in (28) and (29) involve matrices of size $L_{s}\times L_{s}$; since $L_{s}$ is typically $10\sim 30$, the cost of this inversion in marginal. Moreover, if $\rho$ is fixed, the two inverses can be pre-computed [13].

Once again, we tackle (26) with the GMM-denoiser described previously, departing from the standard use of ADMM with a convex regularizer, such as the TV [25] used in [13].

The resulting HS sharpening formulation satisfies Theorem 2. In fact, functions $g_{1}$ and $g_{2}$ (see (22) and (23)) are quadratic, thus closed, proper, and convex. Moreover, using Lemma 2, the denoiser corresponds to the proximity operator of a quadratic function $g_{3}$, which is strongly convex, closed and proper. Finally, since the cost function is the sum of three convex quadratic functions, one of which is strongly convex (defined implicitly by the denoiser), a global minimum is guaranteed to exist. We stress that, although the global objective in (21) is quadratic when $\phi$ is a GMM-based denoiser with fixed weights, and thus fixed ${\bf W}$, we still need to resort to the ADMM to solve the problem because ${\bf W}$ is too large to be handled directly.

6 Results

We compare the results using the GMM-based denoiser from [2], the GMM denoiser herein proposed, and the dictionary-based method from [14] (which, to the best of our knowledge, is the state-of-the-art). We use three different metrics: ERGAS (erreur relative globale adimensionnelle de synthèse), SAM (spectral angle mapper), and SRE (signal to reconstruction-error) [12], in 4 different settings. In the first experiment, we consider sharpening the HS images using the PAN image, at 50dB SNR, on both the HS and PAN images. The second experiment refers to PAN-sharpening as well, with 35dB SNR on the first 43 HS bands and 30dB on the 50 remaining HS bands and the PAN image [14]. Experiments 3 and 4 use the same SNR as experiments 1 and 2, respectively, but the HS bands are sharpened based on 4 MS bands: R, G, B, and near-infrared. In all the experiments, a GMM with 20 components is learned from the PAN image, where the PAN and MS images are generated from the original HS bands, using IKONOS and LANDSAT spectral responses.

Table 1 shows the results on a cropped area of the ROSIS Pavia University and Moffett Field datasets. The three methods have comparable performances, with the proposed one being slightly better. Furthermore, the results supports the hypothesis that the target image has the same spatial structure as the PAN image used to train the GMM, otherwise keeping the posterior mixture weights would not yield good results. Figures 2a–2d show the results of experiment 4 on the Moffett Field data, in terms of visual quality, yet the differences are not noticeable. Another way to visualize the accuracy is to plot the (sorted) errors for each pixel and to notice that the proposed method has a larger number below a given threshold, i.e., in Figure 2e the dashed curve is below the solid one; furthermore, in Figures 2f and 2g, the dashed curve is closer (smaller error) to the dotted one (the true pixel value across all bands) than the solid line.

7 Conclusion

This paper proposed a PnP algorithm, with a GMM-based denoiser, guaranteed to converge to a global minimum of the underlying cost function. The denoiser is scene-adapted and it is the proximity operator of a closed, proper, and convex function; consequently, the standard convergence guarantees of ADMM (in Theorem 1) apply, if the remaining conditions are satisfied. We illustrate the application of the algorithm on a data fusion problem known as HS sharpening.

Experimental results show that the proposed method outperforms another state-of-the-art algorithm based on sparse representations on learned dictionaries [14], for most of the test settings. The proposed scene-adaptation also improves the results over the GMM-based denoiser in [15], giving supporting evidence to the hypothesis that the PAN image constitutes a useful training example, as the HS bands (and latent images) share a similar spatial structure.

Appendix: Proof of Lemma 1

Each ${\bf F}_{i}$ is symmetric, as it is a convex combination of symmetric matrices, which is clear from (16) and because covariances are symmetric. Thus, ${\bf W}$ is also symmetric (see (17)). Consider the eigendecomposition of each {\bf C}_{j}={\bf U}_{j}^{T}\mbox{\boldmath\Sigma}_{j}{\bf U}_{j}, where \mbox{\boldmath\Sigma}_{j}=\mbox{diag}(\varsigma_{1}^{j},...,\varsigma_{n_{p}}^{j}) contains its eigenvalues, in non-increasing order. Then,

[TABLE]

where \mbox{\boldmath\Sigma}_{j}\left(\mbox{\boldmath\Sigma}_{j}+\sigma^{2}\;{\bf I}\right)^{-1} is a diagonal matrix, and thus the eigenvalues of ${\bf C}_{j}({\bf C}_{j}+\sigma^{2}\;{\bf I})^{-1}$ are in

[TABLE]

since $\varsigma_{n_{p}}^{j}>0$ (${\bf C}_{j}$ is positive definite) and $\sigma^{2}>0$.

From (16), each ${\bf F}_{i}$ is a convex combination of matrices, each of which with eigenvalues in $(0,1)$. Weyl’s inequality [26] implies that the eigenvalues of a convex combination of symmetric matrices is bounded below (above) by the same convex combination of the smallest (largest) eigenvalues of those matrices. The eigenvalues of ${\bf F}_{i}$ are thus all in $(0,1)$, i.e., ${\bf F}_{i}$ is positive definite and $\|{\bf F}_{i}\|_{2}<1$.

Finally, following [27], we partition the set of patches into a collection of subsets of non-overlapping patches: $\{\Omega_{j}\subset\{1,..,,N\},\;j=1,...,n_{p}\}$ (the number of subsets of non-overlapping patches is equal to the patch size, due to the assumption of unit stride and periodic boundaries). Using this partition, ${\bf W}$ can be written as

[TABLE]

Since the patches in $\Omega_{j}$ are disjoint, there is a permutation of the image pixels that allows writing ${\bf A}_{j}$ as a block-diagonal matrix, where the blocks are the ${\bf F}_{k}$ matrices, with $k\in\Omega_{j}$. Because the set of eigenvalues of a block-diagonal matrix is the union of the sets of eigenvalues of its blocks, the eigenvalues of each ${\bf A}_{j}$ are bounded similarly as those of the ${\bf F}_{k}$, thus in $(0,1)$. Finally, again using Weyl’s inequality, the eigenvalues of ${\bf W}$ are bounded above (below) by the average of the largest (smallest) eigenvalues of each ${\bf A}_{j}$, thus also in $(0,1)$. $\blacksquare$