Learning Convex Regularizers for Optimal Bayesian Denoising

TL;DR

This paper introduces a data-driven, convex regularizer learning framework for Bayesian denoising that guarantees convergence and adapts easily to different noise levels, achieving near-MMSE performance.

Contribution

It proposes a novel ADMM-based method that learns convex penalty functions via spline parametrization, with proven convergence and noise-adaptive capabilities.

Findings

Achieves denoising performance comparable to MMSE estimators.

Provides a convergent, noise-adaptive learning framework for Bayesian denoising.

Demonstrates effectiveness on sparse and non-sparse Lévy process models.

Abstract

We propose a data-driven algorithm for the maximum a posteriori (MAP) estimation of stochastic processes from noisy observations. The primary statistical properties of the sought signal is specified by the penalty function (i.e., negative logarithm of the prior probability density function). Our alternating direction method of multipliers (ADMM)-based approach translates the estimation task into successive applications of the proximal mapping of the penalty function. Capitalizing on this direct link, we define the proximal operator as a parametric spline curve and optimize the spline coefficients by minimizing the average reconstruction error for a given training set. The key aspects of our learning method are that the associated penalty function is constrained to be convex and the convergence of the ADMM iterations is proven. As a result of these theoretical guarantees, adaptation of the proposed framework to different levels of measurement noise is extremely simple and does not require any retraining. We apply our method to estimation of both sparse and non-sparse models of L\'{e}vy processes for which the minimum mean square error (MMSE) estimators are available. We carry out a single training session and perform comparisons at various signal-to-noise ratio (SNR) values. Simulations illustrate that the performance of our algorithm is practically identical to the one of the MMSE estimator irrespective of the noise power.

Equations

$$\mathrm{dom}T$$

$$\mathrm{ran}T$$

$$\mathrm{gra}T$$

$$T^{-1}\mathbf{u}:=\left\{\mathbf{x}\in \mathcal{H}|\mathbf{u}\in T\mathbf{x}\right\}.$$

$$⟨\mathbf{x}-\mathbf{y},\mathbf{u}-\mathbf{v}⟩\ge 0,\forall (\mathbf{x},\mathbf{u})\in \mathrm{gra}T,\forall (\mathbf{y},\mathbf{v})\in \mathrm{gra}T.$$

$$\sum _{i=1}^n⟨{\mathbf{x}}_{i+1}-{\mathbf{x}}_i,{\mathbf{u}}_i⟩\le 0.$$

$$⟨\mathbf{x}-\mathbf{y},\mathbf{u}-\mathbf{v}⟩\ge \Vert \mathbf{u}-\mathbf{v}{\Vert }^2,\forall (\mathbf{x},\mathbf{u})\in \mathrm{gra}T,\forall (\mathbf{y},\mathbf{v})\in \mathrm{gra}T.$$

$$\partial f(\mathbf{x})=\left\{\mathbf{u}\in \mathcal{H}|⟨\mathbf{y}-\mathbf{x},\mathbf{u}⟩\le f(\mathbf{y})-f(\mathbf{x}),\forall \mathbf{y}\in \mathcal{H}\right\},$$

$${\mathrm{prox}}_f(\mathbf{x})=\underset{\mathbf{u}\in \mathcal{H}}{\mathrm{argmin}}\left\{f(\mathbf{u})+\frac{1}{2}\Vert \mathbf{u}-\mathbf{x}{\Vert }_2^2\right\}.$$

$${\mathrm{prox}}_f={\left(\mathrm{Id}+\partial f\right)}^{-1},$$

$$p_U(\mathbf{u})=\prod _{i=1}^N{p}_U(u_i),$$

$${\Phi }_U(\mathbf{u})=\sum _{i=1}^N{\Phi }_U(u_i),$$

$${\hat{\mathbf{x}}}_{\mathrm{MMSE}}(\mathbf{y})=\mathbf{y}+{\sigma }^2\nabla \log p_Y(\mathbf{y}),$$

$${\hat{\mathbf{x}}}_{\mathrm{MAP}}(\mathbf{y})$$

$$=\underset{\mathbf{x}}{\mathrm{argmax}}\left\{p_{Y|X}(\mathbf{y}|\mathbf{x})p_X(\mathbf{x})\right\}$$

$$=\underset{\mathbf{x}}{\mathrm{argmin}}\left\{\frac{1}{2}\Vert \mathbf{y}-\mathbf{x}{\Vert }_2^2+{\sigma }^2{\Phi }_X(\mathbf{x})\right\}.$$

$$\underset{\mathbf{x}}{\min }\left\{\frac{1}{2}\Vert \mathbf{y}-\mathbf{x}{\Vert }_2^2+{\sigma }^2\sum _{i=1}^N{\Phi }_U([\mathbf{L}\mathbf{x}{]}_i)\right\}.$$

$$\frac{1}{2}\Vert \mathbf{y}-\mathbf{x}{\Vert }_2^2+{\sigma }^2{\Phi }_U(\mathbf{u})-⟨\mathbf{\alpha },\mathbf{L}\mathbf{x}-\mathbf{u}⟩+\frac{\mu }{2}\Vert \mathbf{L}\mathbf{x}-\mathbf{u}{\Vert }_2^2$$

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

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

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

$$T(x)=\sum _{m=-M}^M{c}_m\psi \left(\frac{x}{\Delta }-m\right),$$

$$J(\mathbf{c})=\frac{1}{2}\sum _{\ell =1}^L{\Vert {\mathbf{x}}^{(K)}(\mathbf{c},{\mathbf{y}}_{\ell })-{\mathbf{x}}_{\ell }\Vert }_2^2,$$

$$J(\mathbf{c})=\frac{1}{2}{\Vert {\mathbf{x}}^{(K)}(\mathbf{c},\mathbf{y})-\mathbf{x}\Vert }_2^2.$$

$${\mathbf{c}}^{(i+1)}={\mathbf{c}}^{(i)}-\gamma \nabla J\left({\mathbf{c}}^{(i)}\right).$$

$$\nabla y(\mathbf{x})=\frac{\partial y}{\partial \mathbf{x}}={\left[\begin{array}{cccc}\frac{\partial y}{\partial x_1}& \frac{\partial y}{\partial x_2}& \cdots & \frac{\partial y}{\partial x_m}\end{array}\right]}^T,$$

$$\frac{\partial \mathbf{y}}{\partial \mathbf{x}}=\left[\begin{array}{cccc}\frac{\partial y_1}{\partial \mathbf{x}}& \frac{\partial y_2}{\partial \mathbf{x}}& \cdots & \frac{\partial y_n}{\partial \mathbf{x}}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial y_1}{\partial x_1}& \cdots & \frac{\partial y_n}{\partial x_1}\\ ⋮& \ddots & ⋮\\ \frac{\partial y_1}{\partial x_m}& \cdots & \frac{\partial y_n}{\partial x_m}\end{array}\right].$$

$$\mathbf{M}$$

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

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

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

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

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

$$\nabla J(\mathbf{c})$$

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

Full text

I Introduction

Statistical inference is a central theme to the theory of inverse problems and, in particular, Bayesian methods have been successfully used in several signal processing problems [1, 2, 3]. Among these is the estimation of signals under the additive white Gaussian noise (AWGN) hypothesis, which we shall consider throughout this paper. Conventionally, the unobservable signal is modeled as a random object with a prior probability density function (pdf) and the estimation is performed by assessing the posterior pdf that characterizes the problem statistically. In addition to being important on its own right, this classical problem has recently gained a significant amount of interest. The main reason of the momentum is that Bayesian estimators can be directly integrated—as “denoisers”—into algorithms that are designed for more sophisticated inverse problems [4]. In plain terms, employing a more accurate denoising technique helps one improve the performance of the subsequent reconstruction method. Such ideas have been presented in various applications including deconvolution [5], super-resolved sensing [6], and compressive imaging [7], to name just a few.

The MAP inference is by far the most widely used Bayesian paradigm due its computational convenience [8]. Compatibility of MAP with the regularized least-squares approach is well-understood. For example, by this parallelism, the implicit statistical links between the popular sparsity-based methods [9] and MAP considerations based on generalized Gaussian, Laplace, or hyper-Laplace priors is established [10, 11, 12]. Moreover, recent iterative optimization techniques including (fast) iterative shrinkage/thresholding algorithm ((F)ISTA) [13, 14, 15] and ADMM [16] allow us to handle these type of problems very efficiently.

Fundamentally, the estimation performance of MAP is differentiated by the preferred prior model. When the inherent nature of the underlying signal is (fully or partially) deterministic, identification of the right prior is challenging. Fitting statistical models to such signals (or collections of them) is feasible [17]. Yet, the apparent downside is that the reference pdf, which specifies the inference, can be arbitrary. Even when the signal of interest is purely stochastic and the prior is exactly known, deviations from the initial statistical assumptions is observed [18]. More importantly, mathematical characterization of the MAP estimate (i.e., the maximizer of the posterior pdf) by means of mean-square error (MSE) is available only in limited cases[19]. Hence, algorithms driven by rigorous MAP considerations can still be suboptimal with respect to MSE [20, 21]. These observations necessitate revisiting MAP-like formulations from the perspective of estimation accuracy instead of strict derivations based on the prior model.

I-A Overview of Related Literature

Several works have aimed at improving the performance of MAP. Cho et al. have introduced a nonconvex method to enforce the strict fit between the signal (or its attributes) and their choice of prior distribution. Gribonval has shown that the MMSE can actually be stated as a variational problem that it is in spirit of MAP [22]. Based on the theory of continuous-domain sparse stochastic processes [23], Amini et al. have analyzed the conditions under which the performance of MAP can be MSE-optimal [24]. In [25], Bostan et al. have investigated the algorithmic implications of various prior models for the proximal (or the shrinkage) operator that takes part in the ADMM steps. Accordingly, Kazerouni et al. and Tohidi et al. have demonstrated that MMSE performance can be achieved for certain type of signals if the said proximal operator is replaced with carefully chosen MMSE-type shrinkage functions [26]. Such methods, however, rely on the full knowledge of the prior model, which significantly limits their applicability.

Modification of the proximal operators have also been investigated based on deterministic principles. In particular, motivated by the outstanding success of convolutional neural networks (CNNs) [27], several researchers have used learning-based methods to identify model parameters (thus the proximal). In this regard, Gregor and LeCun [28], and Kamilov and Mansour [29] have considered sparse encoding applications and replaced the soft-thresholding step in (F)ISTA with a learned proximal. Yang et al. have applied learning to the proximal operator of ADMM for improved magnetic resonance (MR) image reconstruction [30]. In [31, 32, 33], learning different shrinkage functions for each iteration is proposed. A variant of these methods is considered by Lefkimmiatis [34]. More relevant to the present context, Samuel and Tappen have learned the model parameters of MAP estimators for continuous-valued Markov random fields (MRFs) [35]. What is common in all these techniques is that the proximal algorithm at hand is trained, which is a nonconvex optimization problem, without any restrictions.

I-B Contributions

We revisit the MAP problem that is cast as the minimization of a quadratic fidelity term regularized by a penalty function. The latter captures the statistics of the original signal. The problem is solved via ADMM by iteratively applying the proximal operator associated with the penalty function. This direct link provides us with the proper framework to formulate and rigorously analyze our method. Our main contributions are summarized as follows:

• $\blacksquare$

Proposal of a new estimator by learning an iteration-independent penalty function that is convex. The convexity constraint is appropriately characterized in terms of the spline coefficients that parameterize the corresponding proximal operator. The learning process optimizes the coefficients so that the mean $\ell_{2}$-normed error between a set of ground-truth signals and the ADMM reconstructions (from their noise-added versions) is minimized.

• $\blacksquare$

Convergence proof of the resultant ADMM scheme based on the above-mentioned convexity confinement. Consequently, the learned penalty function is adjusted from one level of noise to another by a simple scaling operation, eliminating the need for retraining. Furthermore, assuming symmetrically distributed signals, the number of learning parameters is reduced by a half.

• $\blacksquare$

Application of the proposed learning framework on two model signals, namely the Brownian motion and compound Poisson process. The main reason for choosing these models is that their (optimal) MMSE estimations are available for comparison. Furthermore, since these stochastic processes can be decorrelated by the finite difference operator, dictionary learning is no longer needed and we can focus only on the nonlinearity learning. Experiments show that, for a wide range of noise variances, ADMM reconstructions with learned penalty functions are almost identical to the minimum mean square error (MMSE) estimators of these signals. We further demonstrate the practical advantages of the proposed learning scheme over its unconstrained counterpart.

I-C Outline

In the sequel, we provide an overview of the necessary mathematical tools in Section II. In Section III, we present our spline-based parametrization for the proximal operator and formulate the unconstrained version of our algorithm. This is then followed by the introduction of the constraint formulation in terms of the spline coefficients in Section IV. We prove the convergence and the scalability (with respect to noise power) in Section V. Finally, numerical results are illustrated in Section VI where we show that our algorithm achieves the MMSE performance for Lévy processes with different sparsity characteristics.

II Background

II-A Monotone operator theory

We review here some notation and background from convex analysis and monotone operator theory; see [36] for further details. Let us restrict ourselves to the Hilbert space ${\mathcal{H}}={\mathbb{R}}^{d}$, for some dimension $d\geq 1$, equipped with the Euclidean scalar product $\left\langle{\cdot}\,,\,{\cdot}\right\rangle$ and norm $\|\cdot\|_{2}$. The identity operator on ${\mathcal{H}}$ is denoted by $\operatorname{Id}$. Consider a set-valued operator $T:{\mathcal{H}}\rightarrow 2^{{\mathcal{H}}}$ that maps each vector ${\boldsymbol{x}}\in{\mathcal{H}}$ to a set $T{\boldsymbol{x}}\subset{\mathcal{H}}$. The domain, range, and graph of operator $T$ are respectively defined by

[TABLE]

We say that $T$ is single-valued if $T{\boldsymbol{x}}$ has a unique element for all ${\boldsymbol{x}}\in\operatorname{dom}T$. The inverse $T^{-1}$ of $T$ is also a set-valued operator from ${\mathcal{H}}$ to ${2^{{\mathcal{H}}}}$ defined by

[TABLE]

It is straightforward that to see that $\operatorname{dom}T=\operatorname{ran}T^{-1}$ and $\operatorname{ran}T=\operatorname{dom}T^{-1}$. $T$ is called monotone if

[TABLE]

In 1-D, a monotone operator is simply a non-decreasing function. $T$ is maximally monotone if it is monotone and there exists no monotone operator $S$ such that $\operatorname{gra}T\subsetneqq\operatorname{gra}S$. A handy characterization of the maximal monotonicity is given by Minty’s theorem [36, Theorem 21.1].

Theorem 1** (Minty).**

A monotone operator $T:{\mathcal{H}}\rightarrow 2^{{\mathcal{H}}}$ is maximally monotone if and only if $\operatorname{ran}(\operatorname{Id}+T)={\mathcal{H}}$.

For an integer $n\geq 2$, $T:{\mathcal{H}}\rightarrow 2^{{\mathcal{H}}}$ is $n$-cyclically monotone if, for every $n$ points $({\boldsymbol{x}}_{i},{\boldsymbol{u}}_{i})\in\operatorname{gra}T,i=1,\ldots,n$, and for $x_{n+1}=x_{1}$, we have that

[TABLE]

An operator is cyclically monotone if it is $n$-cyclically monotone for all $n\geq 2$. This is a stronger notion of monotonicity because being monotone is equivalent to being $2$-cyclically monotone. Moreover, $T$ is maximally cyclically monotone if it is cyclically monotone and there exists no cyclically monotone operator $S$ such that $\operatorname{gra}T\subsetneqq\operatorname{gra}S$. An operator $T$ is said to be firmly nonexpansive if

[TABLE]

It is not difficult to see that a firmly nonexpansive operator must be both single-valued and monotone.

We denote by $\Gamma_{0}({\mathcal{H}})$ the class of all proper lower-semicontinuous convex functions $f:{\mathcal{H}}\rightarrow(-\infty,+\infty]$. For any proper function $f:{\mathcal{H}}\rightarrow(-\infty,+\infty]$, the subdifferential operator $\partial f:{\mathcal{H}}\rightarrow 2^{{\mathcal{H}}}$ is defined by

[TABLE]

whereas, the proximal operator $\operatorname{prox}_{f}:{\mathcal{H}}\rightarrow 2^{{\mathcal{H}}}$ is given by

[TABLE]

It is remarkable that, when $f\in\Gamma_{0}({\mathcal{H}})$, $\partial f$ is maximally cyclically monotone, $\operatorname{prox}_{f}$ is firmly nonexpansive, and the two operators are related by

[TABLE]

where the right-hand side is also referred to as the resolvent of $\partial f$. Interestingly, any maximally cyclically monotone operator is the subdifferential of some convex function, according to Rockafellar’s theorem [36, Theorem 22.14].

Theorem 2** (Rockafellar).**

$A:{\mathcal{H}}\rightarrow 2^{{\mathcal{H}}}$* is maximally cyclically monotone if and only if there exists $f\in\Gamma_{0}({\mathcal{H}})$ such that $A=\partial f$.*

II-B Denoising problem and ADMM

Let us consider throughout this paper the denoising problem in which a signal ${\boldsymbol{x}}\in{\mathbb{R}}^{N}$ is estimated from its corrupted version ${\boldsymbol{y}}={\boldsymbol{x}}+{\boldsymbol{n}}$, where ${\boldsymbol{n}}$ is assumed to be additive white Gaussian noise (AWGN) of variance $\sigma^{2}$. An estimator of ${\boldsymbol{x}}$ from ${\boldsymbol{y}}$ is denoted by $\hat{{\boldsymbol{x}}}({\boldsymbol{y}})$. We treat ${\boldsymbol{x}}$ as a random vector generated from the joint probability density function (pdf) $p_{X}$. It is assumed that ${\boldsymbol{x}}$ is whitenable by a matrix ${\boldsymbol{L}}\in{\mathbb{R}}^{N\times N}$ such that the transformed vector ${\boldsymbol{u}}={\boldsymbol{L}}{\boldsymbol{x}}$ has identically independently distributed (i.i.d.) entries. The joint pdf $p_{U}$ of the so-called innovation ${\boldsymbol{u}}$ is therefore separable, i.e.,

[TABLE]

where, for convenience, $p_{U}$ is reused to denote the 1-D pdf of each component of ${\boldsymbol{u}}$. We define $\Phi_{U}({\boldsymbol{u}})=-\log p_{U}({\boldsymbol{u}})$ as the penalty function of ${\boldsymbol{u}}$. This function is then separable in the sense that

[TABLE]

where $\Phi_{U}$ is again used to denote the 1-D penalty function of each entry $u_{i}$.

The MMSE estimator, which is optimal if the ultimate goal is to minimize the expected squared error between the estimate $\hat{{\boldsymbol{x}}}$ and the original signal ${\boldsymbol{x}}$, is given by Stein’s formula [19]

[TABLE]

where $p_{Y}$ is the joint pdf of the measurement ${\boldsymbol{y}}$ and $\nabla$ denotes the gradient operator. Despite its elegant expression, the MMSE estimator, in most cases, is computationally intractable since $p_{Y}$ is obtained through a highly dimensional convolution between the prior distribution $p_{X}$ and the Gaussian distribution $g_{\sigma}({\boldsymbol{n}})=(2\pi\sigma^{2})^{-N/2}\exp(-{\left\|{\boldsymbol{n}}\right\|^{2}}/{2\sigma^{2}})$. However, for Lévy processes, which have independent and stationary increments, the MMSE estimator is computable using a message passing algorithm [37].

The maximum a posteriori estimator (MAP) is given by

[TABLE]

where $\Phi_{X}({\boldsymbol{x}})=-\log p_{X}({\boldsymbol{x}})$ is the (nonseparable) penalty function of ${\boldsymbol{x}}$. In other words, the MAP estimator is exactly the proximal operator of $\sigma^{2}\Phi_{X}$. Since ${\boldsymbol{u}}={\boldsymbol{L}}{\boldsymbol{x}}$ is associated with the separable penalty function $\Phi_{U}$, the minimization in (3) can be written as

[TABLE]

This expression of the MAP reconstruction resembles the conventional regularization-based approach in which the transform ${\boldsymbol{L}}$ is designed to sparsify the signal, the penalty function $\Phi_{U}$ is chosen—the typical choice being the $\ell_{1}$-norm—to promote the sparsity of the transform coefficients ${\boldsymbol{u}}$. The parameter $\sigma^{2}$ is set (not necessarily to the noise variance) to trade off the quadratic fidelity term with the regularization term. The optimization problem (4) can be solved efficiently by iterative algorithms such as the alternating direction method of multipliers (ADMM). To that end, we form the augmented Langrangian

[TABLE]

and successively mimimize this functional with respect to each of the variables ${\boldsymbol{x}}$ and ${\boldsymbol{u}}$, while fixing the other one; the Lagrange multiplier ${\boldsymbol{\alpha}}$ is also updated appropriately at each step. In particular, at iteration $k+1$, the updates look like

[TABLE]

Here, ${\boldsymbol{u}}$ and ${\boldsymbol{\alpha}}$ are initialized to be ${\boldsymbol{u}}^{(0)}$ and ${\boldsymbol{\alpha}}^{(0)}$, respectively; ${\boldsymbol{I}}\in{\mathbb{R}}^{N\times N}$ denotes the identity matrix. If the proximal operator $\operatorname{prox}_{\sigma^{2}/\mu\Phi_{U}}$ in (7) is replaced with a general operator $T$, we refer to the above algorithm as the generalized ADMM associated with $T$. When the operator $T$ is separable, i.e., $T({\boldsymbol{u}})=\left(T(u_{1}),\ldots,T(u_{N})\right)$, we refer to the 1-D function $T:{\mathbb{R}}\rightarrow{\mathbb{R}}$ as the shrinkage function; the name comes from the observation that typical proximal operators, such as the soft-thresholding, shrink large values of the input in a pointwise manner to reduce the noise. In what follows, we propose a learning approach to the denoising problem in which the shrinkage function $T$ of the generalized ADMM is optimized in the MMSE sense from data, instead of being engineered as in sparsity-promoting schemes.

III Learning Unconstrained Shrinkage Functions

III-A Learning algorithm

To learn the shrinkage function $T:{\mathbb{R}}\rightarrow{\mathbb{R}}$, we parameterize it via a spline representation:

[TABLE]

where $\psi$ is some kernel (radial basis functions, B-splines, etc.) and $\Delta$ is the sampling step size that defines the distant between consecutive spline knots. We call such function $T$ a shrinkage spline. Consider the generalized ADMM associated with $T$. Fix the transform matrix ${\boldsymbol{L}}$, the penalty parameter $\mu$, and the number of ADMM iterations $K$. The vector ${\boldsymbol{c}}\in{\mathbb{R}}^{2M+1}$ of spline coefficients of $T$ is to be learned by minimizing the following cost function:

[TABLE]

where $\{{\boldsymbol{x}}_{\ell}\}_{\ell=1}^{L}$ is the collection of $L$ ground-truth signals from which the observations $\{{\boldsymbol{y}}_{\ell}\}_{\ell=1}^{L}$ are made and $\left\{{\boldsymbol{x}}^{(K)}({\boldsymbol{y}}_{\ell})\right\}_{\ell=1}^{L}$ are the corresponding reconstructions from $\{{\boldsymbol{y}}_{\ell}\}_{\ell=1}^{L}$ using the generalized ADMM with $K$ iterations. For notational simplicity, from now on we drop the subscript $\ell$ and develop a learning algorithm for a single training example $({\boldsymbol{x}},{\boldsymbol{y}})$ that can be easily generalized to training sets of arbitrary size. The cost function is thus simplified to

[TABLE]

The minimization of this function is performed via a simple gradient descent that is described in Algorithm 1; we call the algorithm unconstrained learning to distinguish it from the constrained learning that will be presented later. As in every neural network, the computation of the gradient of the cost function is performed in a backpropagation manner, which will be detailed in the next section.

III-B Gradient computation

We devise in this section a backpropagation algorithm to evaluate the gradient of the cost function with respect to the spline coefficients of the shrinkage function. We adopt the following convention for matrix calculus: for a function $y:{\mathbb{R}}^{m}\rightarrow{\mathbb{R}}$ of vector variable ${\boldsymbol{x}}$, its gradient is a column vector given by

[TABLE]

whereas, for a vector-valued function ${\boldsymbol{y}}:{\mathbb{R}}^{m}\rightarrow{\mathbb{R}}^{n}$ of vector variable ${\boldsymbol{x}}$, its Jacobian is an $m\times n$ matrix defined by

[TABLE]

We are now ready to compute the gradient of the cost function $J$ with respect to the parameter vector ${\boldsymbol{c}}$. For simplicity, for $k=0,\ldots,K-1$, put

[TABLE]

By using these notations, we concisely write the updates at iteration $k+1$ of the generalized ADMM associated with operator $T$ as

[TABLE]

First, applying the chain rule to (10) yields

[TABLE]

Next, from the updates of the ADMM and by noting that ${\boldsymbol{L}}$ and ${\boldsymbol{y}}$ does not depend on ${\boldsymbol{c}}$, for $k=0,\ldots,K-1$, we get

[TABLE]

[TABLE]

and

[TABLE]

where ${\boldsymbol{D}}^{(k)}=\operatorname{diag}\left(T^{\prime}({\boldsymbol{v}}^{(k)})\right)$ is the diagonal matrix whose entries on the diagonal are the derivatives of $T$ at $\{v^{(k)}_{i}\}_{i=1}^{N}$, and ${\boldsymbol{\Psi}}^{(k)}$ is a matrix defined by

[TABLE]

Proceeding with simple algebraic manipulation, we arrive at

[TABLE]

where

[TABLE]

Finally, by combining (11) with (12) and by noting that ${\partial{\boldsymbol{\alpha}}^{(0)}}/{\partial{\boldsymbol{c}}}={\partial{\boldsymbol{u}}^{(0)}}/{\partial{\boldsymbol{c}}}={\boldsymbol{0}}$, we propose a backpropagation algorithm to compute the gradient of the cost function $J$ with respect to the spline coefficients ${\boldsymbol{c}}$ as described in Algorithm 2. We refer to the generalized ADMM that uses a shrinkage function learned via Algorithm 1 as MMSE-ADMM.

IV Learning Constrained Shrinkage Functions

We propose two constraints for learning the shrinkage functions: firm nonexpansiveness and antisymmetry. The former is motivated by the well-known fact that the proximal operator of a convex function must be firmly nonexpansive [36]; the latter is justified by Theorem 3: symmetrically distributed signals imply antisymmetric proximal operator and vice versa.

Theorem 3**.**

Let $\Phi\in\Gamma_{0}({\mathbb{R}^{N}})$. $\Phi$ is symmetric if and only if $\operatorname{prox}_{\Phi}$ is antisymmetric.

Proof:

See Appendix A. ∎

In order to incorporate the firmly nonexpansive constraint into the learning of spline coefficients, we choose the kernel $\psi$ in the representation (8) to be a B-spline of some integer order. Recall that the B-spline $\beta^{n}$ of integer order $n\geq 0$ is defined recursively as

[TABLE]

These are compactly supported kernels with many good properties [38]. More importantly, as pointed out in Theorem 4, by using B-spline kernels the firm nonexpansiveness of a shrinkage spline is satisfied as long as its coefficients obey a simple linear constraint.

Theorem 4**.**

Let $\Delta>0$ and let $\beta^{n}$ be the B-spline of order $n\geq 1$. If $c$ is a sequence such that $0\leq c_{m}-c_{m-1}\leq\Delta,\forall m\in{\mathbb{Z}}$, then $f=\sum_{m\in{\mathbb{Z}}}c_{m}\beta^{n}(\cdot/\Delta-m)$ is a firmly nonexpansive function.

Proof:

Since $f$ is a 1-D function, it is easy to see that $f$ is firmly expansive if and only if

[TABLE]

We now show (13) by considering 2 different cases.

$n=1$: $\beta^{1}$ is the triangle function, and so $f$ is continuous and piecewise-linear. If $x,y\in[(m-1)\Delta,m\Delta]$ for some $m\in{\mathbb{Z}}$, then

[TABLE]

which implies

[TABLE]

Otherwise, there exist $k,\ell\in{\mathbb{Z}}$ such that $y\in[(k-1)\Delta,k\Delta]$, $x\in[\ell\Delta,(\ell+1)\Delta]$. Then, we write

[TABLE]

By applying (14) to each term of the above sum, we obtain the desired pair of inequalities in (13), which implies the firm nonexpansiveness of $f$.

$n\geq 2$: $\beta^{n}$ is now differentiable and so is $f$. Thus, by using the mean value theorem, (13) is achieved if the derivative $f^{\prime}$ of $f$ is bounded between 0 and 1, which will be shown subsequently. Recall that the derivative of $\beta^{n}$ is equal to the finite difference of $\beta^{n-1}$. In particular,

[TABLE]

Hence, for all $x\in{\mathbb{R}}$,

[TABLE]

Since $0\leq c_{m}-c_{m-1}\leq\Delta,\forall m\in{\mathbb{Z}}$ and since $\beta^{n-1}({x}/{\Delta}+{1}/{2}-m)\geq 0,\forall x\in{\mathbb{R}},m\in{\mathbb{Z}}$, one has the following pair of inequalities for all $x\in{\mathbb{R}}$:

[TABLE]

By using the partition-of-unity property of the B-spline $\beta^{n-1}$, (16) is simplified to

[TABLE]

which finally proves that $f$ is a firmly nonexpansive function. ∎

With the above results, we easily design an algorithm for learning antisymmetric and firmly nonexpansive shrinkage functions. Algorithm 3 is the constrained counterpart of Algorithm 1: the gradient descent is replaced with a projected gradient descent where, at each update, the spline coefficients are projected onto the linear-constraint set described in Theorem 4 (this projection is performed via a quadratic programming). The gradient of the cost function in this case is evaluated through Algorithm 4, which is just slightly modified from Algorithm 2 to adapt to the antisymmetric nature of the shrinkage splines. We refer to the generalized ADMM that uses a shrinkage function learned by Algorithm 3 as MMSE-CADMM, where the letter ‘C’ stands for ‘convex.’ The convexity of this learning scheme will be made clear in Section V.

V Advantages of Adding Constraints

It is clear that imposing the antisymmetric constraint on the shrinkage function reduces the dimension of the optimization problem by a half and therefore substantially reduces the learning time. In this section, we demonstrate, from the theoretical point of view, the two important advantages of imposing the firmly nonexpansive constraint on the shrinkage function: convergence guarantee and scalability with noise level. Thanks to these properties, our learning-based denoiser behaves like a MAP estimator with some convex penalty function that is now different from the conventional penalty function. On the other hand, as experiments later show, the constrained learning scheme nearly achieves the optimal denoising performance of the MMSE estimator.

V-A Convergence guarantee

The following result asserts that the ADMM denoising converges, no matter what the noisy signal is, if it uses a separable firmly nonexpansive operator in the place of the conventional proximal operator. Interestingly, as will be shown in the proof, any separable firmly nonexpansive operator is the proximal operator of a separable convex penalty function. We want to emphasize that the separability is needed to establish this connection, although the reverse statement is known to hold in the multidimensional case [36].

Theorem 5**.**

If $T:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is a 1-D firmly nonexpansive function such that $\operatorname{dom}T={\mathbb{R}}$, then, for every input ${\boldsymbol{y}}\in{\mathbb{R}}^{N}$, the reconstruction sequence $\{{\boldsymbol{x}}^{(k)}\}$ of the generalized ADMM associated with the separable operator $T$ converges to

[TABLE]

as $k\rightarrow\infty$, where $\Phi\in\Gamma_{0}({\mathbb{R}})$ is some 1-D convex function.

Proof:

First, we show that there exists a function $\Phi\in\Gamma({\mathbb{R}})$ such that $T=\operatorname{prox}_{\Phi}$. To that end, let us define

[TABLE]

The firm nonexpansiveness of $T$ then implies

[TABLE]

which means that $S$ is a monotone operator. Furthermore, we have that

[TABLE]

Therefore, $S$ is maximally monotone thanks to Minty’s theorem (Theorem 1). Since $S$ is an operator on ${\mathbb{R}}$, we invoke [36, Thm. 22.18] to deduce that $S$ must also be maximally cyclically monotone. Now, as a consequence of Rockafellar’s theorem (Theorem 2), there exists a function $\Phi\in\Gamma({\mathbb{R}})$ such that $\partial\Phi=S$. By the definition of $S$ and of $\operatorname{prox}_{\Phi}$, we have that

[TABLE]

Next, define the cost function

[TABLE]

Replacing $T$ with $\operatorname{prox}_{\Phi}$, the generalized ADMM associated with $T$ becomes the regular ADMM associated with the above cost function $f$. By using the convexity of the $\ell_{2}$-norm and of the function $\Phi$, it is well known [16, Section 3.2.1] that $f\left({\boldsymbol{x}}^{(k)}\right)\rightarrow p^{*}$ as $k\rightarrow\infty$, where $p^{*}$ is the minimum value of $f$.

Finally, notice that the function $f$ defined in (17) is strongly convex. Thus, there exists a unique minimizer ${\boldsymbol{x}}^{*}\in{\mathbb{R}}^{N}$ such that $f({\boldsymbol{x}}^{*})=p^{*}$. Moreover, ${\boldsymbol{x}}^{*}$ is known to be a strong minimizer [39, Lemma 2.26] in the sense that the convergence of $\left\{f\left({\boldsymbol{x}}^{(k)}\right)\right\}_{k}$ to $f({\boldsymbol{x}}^{*})$ implies the convergence of $\left\{{\boldsymbol{x}}^{(k)}\right\}_{k}$ to ${\boldsymbol{x}}^{*}$. This completes the proof. ∎

V-B Scalability with noise level

In all existing learning schemes, the shrinkage function is learned for a particular level of noise and then applied in the reconstruction of testing signals corrupted with the same level of noise. If the noise variance changes, the shrinkage function must be relearned from scratch, which will create a computational burden on top of the ADMM reconstruction. This drawback is due to the unconstrained learning strategy in which the shrinkage function is not necessarily the proximal operator of any function. Despite its flexibility, arbitrary nonlinearity no longer goes hand-in-hand with a regularization-based minimization (MAP-like estimation). By contrast, our constrained learning scheme maintains a connection with the underlying minimization regularized by a convex penalty function. This strategy allows us to easily adjust the learned shrinkage function from one level of noise to another by simply scaling the corresponding penalty function by the ratio between noise variances like the conventional MAP estimator. Proposition 1 provides a useful formula to extrapolate the proximal operator of a scaled convex function from the proximal operator of the original function.

Proposition 1**.**

For all $f\in\Gamma_{0}({\mathbb{R}^{N}})$,

[TABLE]

Proof:

Recall a basic result in convex analysis [40] that $\partial(\lambda f)=\lambda\partial f$, for all $f\in\Gamma_{0}({\mathbb{R}^{N}})$ and for all $\lambda\geq 0$ . Also recall that the proximal operator of a convex function is the resolvent of the subdifferential operator. Therefore,

[TABLE]

completing the proof. ∎

The next result establishes that all members of the family generated by (18) are firmly nonexpansive when the generator $\operatorname{prox}_{f}$ is replaced with a general firmly nonexpansive operator. It is noteworthy that the result holds in the multidimensional case where a firmly nonexpansive operator is not necessarily the proximal operator of a convex function.

Theorem 6**.**

If $T:{\mathbb{R}^{N}}\rightarrow{\mathbb{R}^{N}}$ is a firmly nonexpansive operator such that $\operatorname{dom}T={\mathbb{R}^{N}}$, then, for all $\lambda\geq 0$, $T_{\lambda}=(\lambda T^{-1}+(1-\lambda)\operatorname{Id})^{-1}$ is firmly nonexpansive and $\operatorname{dom}T_{\lambda}={\mathbb{R}^{N}}$ as well.

Proof:

The claim is trivial for $\lambda=0$. Assume from now on that $\lambda>0$. We first show that $\operatorname{dom}T_{\lambda}={\mathbb{R}^{N}}$. By a straightforward extension of the argument in the proof of Theorem 5 to the multidimensional case, we easily have that the operator $S=T^{-1}-\operatorname{Id}$ is maximally monotone. It follows that $\lambda S$ is also maximally monotone for all $\lambda>0$. By applying Minty’s theorem to the operator $\lambda S$, we obtain

[TABLE]

Next, we show that $T_{\lambda}$ is firmly nonexpansive. Let ${\boldsymbol{x}},{\boldsymbol{y}}\in{\mathbb{R}^{N}}$ and ${\boldsymbol{u}}\in T_{\lambda}({\boldsymbol{x}}),{\boldsymbol{v}}\in T_{\lambda}({\boldsymbol{y}})$. By the definition of $T_{\lambda}$, one readily verifies that

[TABLE]

The firm nonexpansiveness of $T$ yields

[TABLE]

which translates to

[TABLE]

This confirms that $T_{\lambda}$ is a firmly nonexpansive operator. ∎

VI Experimental Results

In this section, we report the experimental denoising results of the two proposed learning schemes: ADMM with unconstrained shrinkage functions learned via Algorithm 1 (denoted MMSE-ADMM) and ADMM with constrained shrinkage functions learned via Algorithm 3 (denoted MMSE-CADMM). Throughout this section, the transform ${\boldsymbol{L}}$ is fixed to be the finite difference operator: $[{\boldsymbol{L}}{\boldsymbol{x}}]_{i}=x_{i}-x_{i-1},\forall i$; the signal length is fixed to $N=100$. Experiments were implemented in MATLAB on the two following types of Lévy processes:

Brownian motion: entries of the increment vector ${\boldsymbol{u}}={\boldsymbol{L}}{\boldsymbol{x}}$ are i.i.d. Gaussian with unit variance: $p_{U}(u)={e^{-{u^{2}}/{2}}}/{\sqrt{2\pi}}$. 2. 2.

Compound Poisson: entries of the increment vector ${\boldsymbol{u}}={\boldsymbol{L}}{\boldsymbol{x}}$ are i.i.d. Bernoulli-Gaussian: $p_{U}(u)=(1-e^{-\lambda})\,{e^{-{u^{2}}/{2}}}/{\sqrt{2\pi}}+e^{-\lambda}\delta(u)$, where $\delta$ is the Dirac impulse and $\lambda=0.6$ is fixed. This results in a piecewise-constant signal ${\boldsymbol{x}}$ with Gaussian jumps.

Specific realizations of these processes and their corrupted versions with AWGN of variance $\sigma^{2}=1$ are shown in Fig. 1.

VI-A Denoising performance

The same parameters were chosen for both learning schemes (constrained and unconstrained). In particular, for each type of processes, a set of 500 signals was used for training and another set of 500 signals was used for testing. The number of ADMM layers (iterations) was set to $K=10$; the penalty parameter of the augmented Lagrangian was set to $\mu=2$. The shrinkage function was represented with the cubic B-spline:

[TABLE]

The spline coefficients $\{c_{m}\}$ were located uniformly in the dynamic range of ${\boldsymbol{u}}={\boldsymbol{L}}{\boldsymbol{x}}$ with sampling step $\Delta=\sigma/2$, which is dependent on the noise level. Learning was performed by running either Algorithm 1 or Algorithm 3 for 1000 iterations with learning rate $\gamma=2\times 10^{-4}$. The shrinkage function was always initialized with the identity line, which corresponds to $c^{(0)}_{m}=m$ for all $m$.

The denoising performances were numerically evaluated by the signal-to-noise ratio (SNR) improvement that is defined by $\Delta{\rm SNR}\,[\rm dB]=10\log_{10}\left({\|\hat{{\boldsymbol{x}}}-{\boldsymbol{x}}\|^{2}_{2}}/{\|{\boldsymbol{y}}-{\boldsymbol{x}}\|^{2}_{2}}\right).$ We compare the results of MMSE-ADMM and MMSE-CADMM against the following reconstruction methods:

MMSE: This is the optimal estimator (in the MSE sense) and is obtained through a message-passing algorithm [37]. 2. 2.

LMMSE (Linear MMSE): The best linear estimation is obtained by applying the Wiener filter to the noisy observation: $\hat{{\boldsymbol{x}}}_{\rm LMMSE}=({\boldsymbol{I}}+\sigma^{2}{\boldsymbol{L}}^{T}{\boldsymbol{L}})^{-1}{\boldsymbol{y}}$. This is also the least-square solution with $\ell_{2}$ (Tikhonov-like) regularization. 3. 3.

TV (Total Variation) [9]: This estimator is obtained with an $\ell_{1}$ regularizer whose proximal operator is simply a soft-thresholding: $T_{\lambda}(u)=1_{\{|u|>\lambda\}}{\rm sign}(u)(|u|-\lambda)$. In our experiments, the regularization parameter $\lambda$ is optimized for each signal.

The denoising performances of MMSE-ADMM and MMSE-CADMM for various noise variances between $10^{-{1}/{2}}$ and $10^{{1}/{2}}$ are reported in Figs. 2 and 3, respectively. It is remarkable that both MMSE-ADMM and MMSE-CADMM curves are almost identical to the optimal MMSE curve (with the largest gap being about 0.1 dB) and significantly outperform TV, for both types of signals, and LMMSE, for compound Poisson processes. Note that, for Brownian motions, LMMSE and MMSE are the same. The unconstrained and constrained shrinkage functions learned for three different levels of noise are illustrated in Figs. 4 and 5, respectively. As can be seen in Fig. 4, the unconstrained learning might result in non-monotonic curves, which cannot be the proximal operators of any penalty functions, according to [31, Proposition 1]. By contrast, the antisymmetric and firmly nonexpansive curves in Fig. 5 are the proximal operators of the symmetric and convex penalty functions that are plotted in Fig. 6. These functions were numerically obtained by integrating $\partial\Phi=(T-\operatorname{Id})^{-1}$, where $T$ is the learned shrinkage function.

VI-B Constrained versus unconstrained learning

To demonstrate the benefits of the constrained learning over the unconstrained one, we compare their denoising performances for 9 different levels of noise as before, but this time only the shrinkage function $T$ for $\sigma^{2}=1$ was learned. For constrained learning, the shrinkage function with respect to another noise variance $\sigma^{2}$ was numerically computed by using the formula

[TABLE]

For unconstrained learning, these computations are prohibited, and so the learned shrinkage function for $\sigma^{2}=1$ was used for all the other noise levels. The results were illustrated in Fig. 7. It is noticeable that MMSE-CADMM is much better than MMSE-ADMM and, surprisingly, almost as good as the optimal MMSE for all levels of noise, even though the (constrained) learning was performed only once. In other words, the experiments suggest that the proposed MMSE-CADMM combines desired properties of the MAP and MMSE estimators: fast implementation and scalability with noise variance of MAP and optimality of MMSE.

Another advantage of the constrained learning is its convergence guarantee that is associated with the minimization of an underlying cost function (as mentioned in Theorem 5), which does not necessarily exist in the case of unconstrained learning. Figs. 8 and 9 illustrates the reconstructions of a Brownian motion and a compound Poisson signal, respectively, from their noisy measurements using MMSE-CADMM, and the convergences of the corresponding cost functions. Experiments also show that the constrained learning is much more stable to the number of ADMM iterations used in the testing phase ($K_{\rm test}$) when it is different from the number of ADMM iterations used in the training phase ($K_{\rm train}$). Fig. 10 demonstrates this observation by plotting the average SNRs of denoising compound Poisson signals using MMSE-ADMM and MMSE-CADMM against $K_{\rm test}$ ranging from 2 to 50. In this experiment, both constrained and unconstrained learnings were performed with $K_{\rm train}=2$ and $\sigma^{2}=10$. It can be seen from the plot that, when $K_{\rm test}$ increases, the SNR of MMSE-ADMM tends to decrease and fluctuate significantly, while the SNR of MMSE-CADMM tends to improve and converge.

VII Conclusion

We have developed in this paper a learning scheme for signal denoising using ADMM in which a single (iteration-independent) shrinkage function is constrained to be antisymmetric firmly-nonexpansive and learned from data via a simple projected gradient descent to minimize the reconstruction error. This constrained shrinkage function is proved to be the proximal operator of a symmetric convex penalty function. Imposing constraints on the shrinkage function gains several striking advantages: the antisymmetry reduces the number of learning parameters by a half, while the firm nonexpansiveness guarantees the convergence of ADMM, as well as the scalability with noise level. Yet, the denoising performance of the proposed learning scheme is empirically identical to the optimal MMSE estimators for the two types of Lévy processes in a wide range of noise variances. Our experiments also demonstrate that learning the convex penalty function for one level of noise (via learning its proximal operator) and then scaling it for other noise levels yields equivalent performances to those of direct leaning for all noise levels. This property opens up an opportunity to vastly improve the robustness and generalization ability of learning schemes. Potential directions for future research include extension of the proposed framework to general inverse problems of the form ${\boldsymbol{y}}={\boldsymbol{H}}{\boldsymbol{x}}+{\boldsymbol{n}}$ as well as to multidimensional signals. Another issue worth investigating is the joint learning of the shrinkage function and the decorrelation (sparsifying) transform ${\boldsymbol{L}}$ from real data, like images, whose statistics are unknown.

Appendix A Proof of Theorem 3

We first recall that

[TABLE]

Assume for now that $\Phi$ is symmetric. Fix ${\boldsymbol{x}}\in{\mathbb{R}^{N}}$ and let ${\boldsymbol{u}}=\operatorname{prox}_{\Phi}({\boldsymbol{x}})$, ${\boldsymbol{v}}=\operatorname{prox}_{\Phi}(-{\boldsymbol{x}})$. We need to show that ${\boldsymbol{u}}=-{\boldsymbol{v}}$. From (19), we have that

[TABLE]

By the definition of the subdifferential operator, we obtain the following inequalities:

[TABLE]

which, by the symmetry of $\Phi$, can be further simplified to

[TABLE]

Adding these inequalities yields $\left\|{\boldsymbol{u}}+{\boldsymbol{v}}\right\|^{2}_{2}\leq 0$, or ${\boldsymbol{u}}=-{\boldsymbol{v}}$.

Assume conversely that $\operatorname{prox}_{\Phi}$ is antisymmetric. We first show that ${\boldsymbol{u}}\in\partial\Phi({\boldsymbol{x}})$ is equivalent to $-{\boldsymbol{u}}\in\partial\Phi(-{\boldsymbol{x}})$. Indeed, by using (19) and from the antisymmetry of $\operatorname{prox}_{\Phi}$,

[TABLE]

Furthermore, $\operatorname{prox}_{\Phi}({\boldsymbol{0}})={\boldsymbol{0}}$ due to the antisymmetry. Since $\partial\Phi({\boldsymbol{0}})=\operatorname{prox}^{-1}_{\Phi}({\boldsymbol{0}})$, it must be that ${\boldsymbol{0}}\in\partial\Phi({\boldsymbol{0}})$. Let $G=\operatorname{gra}(\partial\Phi)$ and choose $({\boldsymbol{x}}_{0},{\boldsymbol{u}}_{0})=({\boldsymbol{0}},{\boldsymbol{0}})\in G$. Consider the Rockafellar anti-derivative [40] of $\partial\Phi$:

[TABLE]

It is well known [36, Proposition 22.15] that $f\in\Gamma_{0}({\mathbb{R}^{N}})$ and $\partial f=\partial\Phi$. Therefore, we can invoke [36, Proposition 22.15] to deduce that $\Phi=f+c$, for some constant $c\in{\mathbb{R}}$. To show the symmetry of $\Phi$, it suffices to show the symmetry of $f$. From (24) and by the symmetry of $G$, $f(-{\boldsymbol{x}})$ is equal to

[TABLE]

which shows that $f$ is symmetric, completing the proof.