Charles-Alban Deledalle
email: charles-alban.deledalle@math.u-bordeaux.fr

Shibin Parameswaran

Truong Q Nguyen

Image denoising with generalized Gaussian mixture model patch priors

Keywords: Generalized Gaussian distribution, Mixture models, Image denoising, Patch priors AMS subject classifications. 68U10, 62H35, 94A08

Patch priors have become an important component of image restoration. A powerful approach in this category of restoration algorithms is the popular Expected Patch Log-Likelihood (EPLL) algorithm. EPLL uses a Gaussian mixture model (GMM) prior learned on clean image patches as a way to regularize degraded patches. In this paper, we show that a generalized Gaussian mixture model (GGMM) captures the underlying distribution of patches better than a GMM. Even though GGMM is a powerful prior to combine with EPLL, the non-Gaussianity of its components presents major challenges to be applied to a computationally intensive process of image restoration. Specifically, each patch has to undergo a patch classification step and a shrinkage step. These two steps can be efficiently solved with a GMM prior but are computationally impractical when using a GGMM prior. In this paper, we provide approximations and computational recipes for fast evaluation of these two steps, so that EPLL can embed a GGMM prior on an image with more than tens of thousands of patches. Our main contribution is to analyze the accuracy of our approximations based on thorough theoretical analysis. Our evaluations indicate that the GGMM prior is consistently a better fit for modeling image patch distribution and performs better on average in image denoising task.

1. Introduction. Image restoration is the process of recovering the underlying clean image from its degraded or corrupted observation(s). The images captured by common imaging systems often contain corruptions such as noise, optical or motion blur due to sensor limitations and/or environmental conditions. For this reason, image restoration algorithms have widespread applications in medical imaging, satellite imaging, surveillance, and general consumer imaging applications. Priors on natural images play an important role in image restoration algorithms. Image priors are used to denoise or regularize ill-posed restoration problems such as deblurring and super-resolution, to name just a few. Early attempts in designing image priors relied on modeling local pixel gradients with Gibbs distributions [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF], Laplacian distributions (total variation) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Yang | An efficient tvl1 algorithm for deblurring multichannel images corrupted by impulsive noise[END_REF], hyper-Laplacian distribution [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF], generalized Gaussian distribution [START_REF] Bouman | A generalized gaussian image model for edge-preserving map estimation[END_REF], or Gaussian mixture models [START_REF] Fergus | Removing camera shake from a single photograph[END_REF]. Concurrently, priors have also been designed by modeling coefficients of an image in a transformed domain using generalized Gaussian [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF][START_REF] Do | Wavelet-based texture retrieval using generalized gaussian density and kullback-leibler distance[END_REF] or scaled mixture of Gaussian [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF] priors for wavelet or curvelet coefficients [START_REF] Boubchir | Multivariate statistical modeling of images with the curvelet transform[END_REF]. Alternatively, modeling the distribution of patches of an image (i.e., small windows usually of size 8×8) has proven to be a powerful solution. In particular, popular patch techniques rely on non-local self-similarity [START_REF] Buades | A non-local algorithm for image denoising[END_REF], fields of experts [START_REF] Roth | Fields of experts: A framework for learning image priors[END_REF], learned patch dictionaries [START_REF] Aharon | k-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF][START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF], sparse or low-rank properties of stacks of similar patches [START_REF] Dabov | Image denoising by sparse 3-d transformdomain collaborative filtering[END_REF][START_REF] Deledalle | Image denoising with patch based pca: local versus global[END_REF][START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF][START_REF] Ji | Robust video restoration by joint sparse and low rank matrix approximation[END_REF], patch re-occurrence priors [START_REF] Michaeli | Blind deblurring using internal patch recurrence[END_REF], or more recently mixture models patch priors [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF][START_REF] Van Den | The student-t mixture as a natural image patch prior with application to image compression[END_REF][START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising (HDMI)[END_REF][START_REF] Singh | Adapted non-gaussian mixture model for image denoising[END_REF][START_REF] Niknejad | Class-specific image denoising using importance sampling[END_REF].

Of these approaches, a successful approach introduced by Zoran and Weiss [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] is to model patches of clean natural images using Gaussian Mixture Model (GMM) priors. The agility of this model lies in the fact that a prior learned on clean image patches can be effectively employed to restore a wide range of inverse problems. It is also easily extendable to include other constraints such as sparsity or multi-resolution patches [START_REF] Sulam | Expected patch log likelihood with a sparse prior[END_REF][START_REF] Papyan | Multi-scale patch-based image restoration[END_REF]. The use of GMMs for patch priors make these methods computationally tractable and flexible. Although GMM patch prior is effective and popular, in this article, we argue that a generalized Gaussian mixture model (GGMM) is a better fit for image patch prior modeling. Compared to a Gaussian model, a generalized Gaussian distribution (GGD) has an extra degree of freedom controlling the shape of the distribution and it encompasses Gaussian and Laplacian models.

Beyond image restoration tasks, GGDs have been used in several different fields of image and signal processing, including watermark detection [START_REF] Cheng | Robust optimum detection of transform domain multiplicative watermarks[END_REF], texture retrieval [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian density and kullback-leibler distance[END_REF], voice activity detection [START_REF] Gazor | Speech probability distribution[END_REF] and MP3 audio encoding [START_REF] Dominguez-Molina | A practical procedure to estimate the shape parameter in the generalized gaussian distribution[END_REF], to cite just a few. In these tasks, GGDs are used to characterize or model the prior distribution of clean signals, for instance, from their DCT coefficients or frequency subbands for natural images [START_REF] Westerink | Subband coding of color images[END_REF][START_REF] Tanabe | Subband image coding using entropy-coded quantization over noisy channels[END_REF][START_REF] Müller | Distribution shape of two-dimensional dct coefficients of natural images[END_REF][START_REF] Aiazzi | Estimation based on entropy matching for generalized gaussian pdf modeling[END_REF] or videos [START_REF] Sharifi | Estimation of shape parameter for generalized gaussian distributions in subband decompositions of video[END_REF], gradients for X-ray images [START_REF] Bouman | A generalized gaussian image model for edge-preserving map estimation[END_REF], wavelet coefficients for natural [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Cheng | Robust optimum detection of transform domain multiplicative watermarks[END_REF], textured [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian density and kullback-leibler distance[END_REF], or ultrasound images [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF], tangential wavelet coefficients for three-dimensional mesh data [START_REF] Lavu | Estimation-quantization geometry coding using normal meshes[END_REF], short time windows for speech signals [START_REF] Kokkinakis | Exponent parameter estimation for generalized gaussian probability density functions with application to speech modeling[END_REF] or frequency subbands for speech [START_REF] Gazor | Speech probability distribution[END_REF] or audio signals [START_REF] Dominguez-Molina | A practical procedure to estimate the shape parameter in the generalized gaussian distribution[END_REF].

In this paper, we go one step further and use multi-variate GGD with one scale and one shape parameter per dimension. The superior patch prior modeling capability of such a GGMM over a GMM is illustrated in Figure 1. The figure shows histograms of six orthogonal 1-D projections of subset of clean patches onto the eigenvectors of the covariance matrix of a single component of the GMM. To illustrate the difference in the shapes (ν) and scales (λ) of the distributions of each dimension, we have chosen a few projections corresponding to both the most and the least significant eigenvalues. It can be seen that GGD is a better fit on the obtained histograms than a Gaussian model. Additionally, different dimensions of the patch follow a different GGD (i.e., has a different shape and scale parameter). Hence, it does not suffice to model all the feature dimensions of a given cluster of patches as Laplacian or Gaussian. Therefore, we propose to model patch priors as GGMM distributed with a separate shape and scale parameters for each feature dimension of a GGD component. This differs from the recent related approach in [START_REF] Niknejad | Class-specific image denoising using importance sampling[END_REF] that considered GGMM where each component has a fixed shape parameter for all directions.

Contributions. The goal of this paper is to measure the improvements obtained in image denoising tasks by incorporating a GGMM in EPLL algorithm. Unlike [START_REF] Niknejad | Class-specific image denoising using importance sampling[END_REF], that incorporates a GGMM prior in a posterior mean estimator based on importance sampling, we directly extend the maximum a posteriori formulation of Zoran and Weiss [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] for the case of GGMM priors. While such a GGMM prior has the ability to capture the underlying distribution of clean patches more closely, we will show that it introduces two major computational challenges in this case. The first one can be thought of as a classification task in which a noisy patch is assigned to one of the components of the mixture. The second one corresponds to an estimation task where a noisy patch is denoised given that it belongs to one of the components of the mixture. Due to the interaction of the noise distribution with the GGD prior, we first show that these two tasks lead to a group of one-dimensional integration and optimization problems, Figure 1. Histograms of the projection of 200, 000 clean patches on 6 eigenvectors j = 1, 2, 3, 4, 62 and 63 of the covariance matrix of one component k of the mixture (with weight w k = 1.3%). The contribution of each clean patch in the histograms is given by its membership values onto this component k (as obtained during the E-Step of EM). For each histogram, a generalized Gaussian distribution was adjusted by estimating the parameters λ and ν by moment estimation (as obtained during the M-Step of our modified EM exposed in Subsection 3.1). For comparisons, we have also provided illustrations of the best fit obtained with a Gaussian distribution.

respectively. Specifically, for x ∈ R, these problems are of the following forms

R exp - (t -x) 2 2σ 2 - |t| ν λ ν ν dt and argmin t∈R (t -x) 2 2σ 2 + |t| ν λ ν ν , (1)
for some ν > 0, σ > 0 and λ ν > 0. In general, they do not admit closed-form solutions but some particular solutions or approximations have been derived for the estimation/optimization problem [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Chaux | A variational formulation for framebased inverse problems[END_REF]. By contrast, up to our knowledge, little is known for approximating the classification/integration one (only crude approximations were proposed in [START_REF] Soury | New results on the sum of two generalized gaussian random variables[END_REF]).

Our contributions are both theoretical-and application-oriented. The major contribution of this paper, which is and theoretical in nature, is to develop an accurate approximation for the classification/integration problem. In particular, we show that our approximation error vanishes for x → 0 and x → ±∞ when ν = 1, see Theorem 1 and Theorem 2. We next generalize this result for 2 3 < ν < 2 in Theorem 3 and Theorem 4. On top of that, we prove that the two problems enjoy some important desired properties in Proposition 2 and Proposition 3. These theoretical results allow the two quantities to be approximated by functions that can be quickly evaluated in order to be incorporated in fast algorithms. Our last contribution is experimental and concerns the performance evaluation of the proposed model in image denoising scenario. For reproducibility, we have released our implementation at https://bitbucket.org/cdeledalle/ggmm-epll.

Potential impacts beyond image denoising. It is important to note that the two main contributions presented in this work, namely approximations for classification and estimation problems, are general techniques that are relevant to a wider area of research problems than image restoration. In particular, our contributions apply to any problems where (i) the underlying clean data are modeled by a GGD or a GGMM, whereas (ii) the observed samples are corrupted by Gaussian noise. They are especially relevant in machine learning scenarios where a GGMM is trained on clean data but data provided during testing time is noisy. That is, our approximations can be used to extend the applicability of the aforementioned clean GGD/GGM based approaches to the less than ideal testing scenario where the data is corrupted by noise. For instance, using the techniques introduced in this paper, one could directly use the GGD based voice activity model of [START_REF] Gazor | Speech probability distribution[END_REF] into the likelihood ratio test based detector of [START_REF] Gazor | A soft voice activity detector based on a laplacian-gaussian model[END_REF] (which relies on solutions of the integration problem but was limited to Laplacian distributions). In fact, we suspect that many studies in signal processing may have limited themselves to Gaussian or Laplacian signal priors because of the complicated integration problem arising from the intricate interaction of GGDs with Gaussian noise. In this paper, we demonstrate that this difficulty can be efficiently overcome with our approximations. This leads us to believe that the impact of the approaches presented in this paper will not only be useful for image restoration but also aid a wider field of general signal processing applications.

Organization. After explaining the considered patch prior based restoration framework in Section 2, we derive our GGMM based restoration scheme in Section 3. The approximations of the classification and estimation problems are studied in Section 4 and Section 5, respectively. Finally, we present numerical experiments and results in Section 6.

Background.

In this section we provide a detailed overview of the use of patch-based priors in Expected Patch Log-Likelihood (EPLL) framework and its usage under GMM priors.

2.1. Image restoration with patch based priors. We consider the problem of estimating an image u ∈ R N (N is the number of pixels) from noisy linear observations v = Au + w, where A : R N → R M is a linear operator and w ∈ R M is a noise component assumed to be white and Gaussian with variance σ 2 . In this paper, we will focus on standard denoising problems where A is the identity matrix, but in more general settings, it can account for loss of information such as blurring. Typical examples for operator A are: a low pass filter (for deconvolution), a masking operator (for inpainting), or a projection on a random subspace (for compressive sensing). To reduce noise and stabilize the inversion of A, some prior information is used for the estimation of u. Recent techniques [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Sulam | Expected patch log likelihood with a sparse prior[END_REF] include this prior information as a model for the distribution of patches found in natural clean images. We consider the EPLL framework [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] that restores an image by maximum a posteriori estimation over all patches, corresponding to the following minimization problem: argmin

u∈R n P 2σ 2 ||Au -v|| 2 - N i=1 log p (P i u) (2)
where P i : R N → R P is the linear operator extracting a patch with P pixels centered at the pixel with location i (typically, P = 8×8), and p(.) is the a priori probability density function (i.e., the statistical model of noiseless patches in natural images). Since i scans all of the N pixels of the image, all patches contribute to the loss and many patches overlap. Allowing for overlapping is important because otherwise there would appear blocking artifacts. While the first term in eq. ( 2) ensures that Au is close to the observations v (this term is the negative log-likelihood under the white Gaussian noise assumption), the second term regularizes the solution u by favoring an image such that all of its patches fit the prior model of patches in natural images.

Optimization with half-quadratic splitting. Problem ( 2) is a large optimization problem where A couples all unknown pixel values of u and the patch prior is often chosen nonconvex. Our method follow the choice made by EPLL of using a classical technique, known as half-quadratic splitting [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF][START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF], that introduces N auxiliary unknown vectors z i ∈ R P , and alternatively consider the penalized optimization problem, for β > 0, as argmin

u∈R n z 1 ,...,z N ∈R P P 2σ 2 ||Au -v|| 2 + β 2 i∈I ||P i u -z i || 2 - i∈I log p (z i ) . (3)
When β → ∞, the problem (3) is equivalent to the original problem [START_REF] Aharon | k-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. In practice, an increasing sequence of β is considered, and the optimization is performed by alternating between the minimization for u and z i . Though little is known about the convergence of this algorithm, few iterations produce remarkable results, in practice. We follow the EPLL settings prescribed in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] by performing 5 iterations of this algorithm with parameter β set to 1 σ 2 {1, 4, 8, 16, 32} for each iteration, respectively. The algorithm is initialized using û = v for the first estimate.

Minimization with respect to u. Considering all z i to be fixed, optimizing (3) for u corresponds to solving a linear inverse problem with a Tikhonov regularization. It has an explicit solution known as the linear minimum mean square estimator (or often referred to as Wiener filtering) which is obtained as: û = argmin

u∈R n P 2σ 2 ||Au -v|| 2 + β 2 i∈I ||P i u -ẑi || 2 = A t A + βσ 2 P i∈I P t i P i -1 A t v + βσ 2 P i∈I P t i ẑi , (4)
where P t i P i is a diagonal matrix whose i-th diagonal element corresponds to the number of patches overlapping the pixel of index i.

Minimization with respect to z i . Considering u to be fixed, optimizing (3) for z i leads to: ẑi ← argmin

z i ∈R P β 2 ||z i -z i || 2 -log p(z i ) where zi = P i û , (5)
which corresponds to the maximum a posterior (MAP) denoising problem under the patch prior p of a patch zi contaminated by Gaussian noise with variance 1/β. The solution of this optimization problem strongly depends on the properties of the chosen patch prior.

Our algorithm will follow the exact same procedure as EPLL by alternating between eq. (4) 1 and (5). In our proposed method, we will be using a generalized Gaussian mixture model (GGMM) to represent patch prior. The general scheme we will adopt to solve eq. ( 5) under GGMM prior is inspired from the one proposed by Zoran & Weiss [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] in the simpler case of Gaussian mixture model (GMM) prior. For this reason, we will now introduce this simpler case of GMM prior before exposing the technical challenges arising from the use of GGMM prior in Section 3.

2.2. Patch denoising with GMM priors. The authors of [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF] suggested using a zeromean Gaussian mixture model (GMM) prior 2 , that, for any patch z ∈ R P , is given by

p (z) = K k=1 w k N P (z; 0 P , Σ k ) (6)
where K is the number of components, w k > 0 are weights such that k w k = 1, and N P (0 P , Σ k ) denotes the multi-variate Gaussian distribution with zero-mean and covariance Σ k ∈ R P ×P . A K-component GMM prior models image patches as being spread over K clusters that have ellipsoid shapes where each coefficient (of each component) follows a Gaussian distribution, i.e., bell-shaped with small tails. In [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], the parameters w k and Σ k of the GMM are learned using the Expectation Maximization algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] on a dataset of 2 million clean patches of size 8 × 8 pixels that are randomly extracted from the training images of the Berkeley Segmentation Database (BSDS) [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. The GMM learned in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] has K = 200 zero-mean Gaussian mixture components.

Due to the multi-modality of the GMM prior, introducing this prior in eq. ( 5) makes the optimization problem highly non-convex:

ẑ ← argmin z∈R P β 2 ||z -z|| 2 -log K k=1 w k N (z; 0 P , Σ k ) . (7)
To circumvent this issue, the EPLL framework introduced by Zoran et al. [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] uses an approximation 3 . In a nutshell, the approximated approach of EPLL provides a way to avoid 1 Since our study focuses only on denoising, we will consider A = IdN . 2 To enforce the zero-mean assumption, patches are first centered on zero, then denoised using eq. ( 5), and, finally, their initial means are added back. In fact, one can show that it corresponds to modeling p(z -z) with a GMM where zj = 1 P i zi for all 1 j P . 3 An alternative investigated in [START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF] is to replace the MAP denoising problem in (5) by the minimum mean square error (MMSE) estimator (a.k.a., posterior mean). The MMSE estimator is defined as an integration problem and has a closed-form solution in case of GMM priors. In our experimental scenarios, this estimator did not lead to significative improvements compared to MAP and we thus did not pursue this idea. the intractability of mixture models, thus making them available to model patch priors [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. Figure 2 provides an illustration of the EPLL framework and the steps involved in solving the optimization problem [START_REF] Bouman | A generalized gaussian image model for edge-preserving map estimation[END_REF] namely:

(i) compute the posterior p(k|z) 4 of each Gaussian component, k = 1, 2 . . . K for the given noisy patch z with assumed noise variance of 1/β, (ii) select the component k that best explains the given patch z, (iii) perform whitening by projecting z over the main directions of that cluster (given by the eigenvectors of Σ k ), and (iv) apply a linear shrinkage on the coefficients with respect to the noise variance 1/β and the spread of the cluster (encoded by the eigenvalues). The details of each of these steps will be discussed in Subsection 3.2 as part of the development of the proposed model which is more general.

In this paper, we suggest using a mixture of generalized Gaussian distributions that will enable image patches to be spread over clusters that are bell shaped in some directions but can be peaky with large tails in others. While the use of GMM priors leads to piece-wise linear estimator (PLE) as a function of z (see [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF]), our GGMM prior will lead to a piecewise non-linear shrinkage estimator.

3. Generalized Gaussian Mixture Models. In this paper, we aim to learn K orthogonal transforms such that each of them can map a subset (cluster) of clean patches into independent zero-mean coefficients. Instead of assuming the coefficient distributions to be bell shaped, we consider that both the scale and the shape of these distributions may vary from one coordinate to another (within the same transform). Our motivation to assume such a highly flexible model is based on the observation illustrated in Figure 1. Given one of such transform and its corresponding cluster of patches, we have displayed the histogram of the patch coefficients for six different coordinates. It can be clearly observed that the shape of the distribution varies depending on the coordinate. Some of them are peaky with heavy tails, and, therefore, would not be faithfully captured by a Gaussian distribution, as done in EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. By contrast, some others have a bell shape, and so would not be captured properly by a peaky and heavy tailed distribution, as done for instance by sparse models [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Aharon | k-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF][START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF][START_REF] Sulam | Expected patch log likelihood with a sparse prior[END_REF]. This shows that one cannot simultaneously decorrelate and sparsify a cluster of clean patches for all coordinates. Since some of the coordinates reveal sparsity while some others reveal Gaussianity, we propose to use a more flexible model that can capture such variations. We propose using a multi-variate zero-mean generalized Gaussian mixture model (GGMM)

p (z) = K k=1 w k G(z; 0 P , Σ k , ν k ) (8)
where K is the number of components and w k > 0 are weights such that k w k = 1. The notation G(0 P , Σ, ν) denotes the P -dimensional generalized Gaussian distribution (GGD) with zero-mean, covariance Σ ∈ R P ×P (symmetric positive definite) and shape parameter ν ∈ R P , whose expression is

G(z; 0 P , Σ, ν) = K 2|Σ ν | 1/2 exp -||Σ -1/2 ν z|| ν ν with ||x|| ν ν = P j=1 |x j | ν j , (9)
where

K = P j=1 ν j Γ(1/ν j )
and

Σ 1/2 ν = Σ 1/2      Γ(1/ν 1 ) Γ(3/ν 1 )
. . .

Γ(1/ν P ) Γ(3/ν P )      . ( 10
)
Denoting the eigen decomposition of matrix Σ by Σ = U ΛU t such that U ∈ R P ×P is unitary and Λ = diag(λ 1 , λ 2 , . . . , λ P ) 2 is diagonal with positive diagonal elements λ 2 j , Σ 1/2 in the above expression is defined as Σ 1/2 = U Λ 1/2 and Σ -1/2 = Λ -1/2 U t is its inverse.

When ν is a constant vector with all entries equal to ν j = 2, G(0 P , Σ, ν) is the multivariate Gaussian distribution N (0 P , Σ) (as used in EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]). When all ν j = 1, it is the multi-variate Laplacian distribution and the subsequent GGMM is a Laplacian Mixture Model (LMM). When all ν j < 1, it is the multi-variate hyper-Laplacian distribution and the subsequent GGMM is a hyper-Laplacian Mixture Model (HLMM). Choosing K = 1 with a constant vector ν corresponds to ν regularization [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF]. But as motivated earlier, unlike classical multivariate GGD models [START_REF] Boubchir | Multivariate statistical modeling of images with the curvelet transform[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF][START_REF] Niknejad | Class-specific image denoising using importance sampling[END_REF], we allow for the entries of ν to vary from one coordinate j to another. To the best of our knowledge, the proposed work is the first one to consider this fully flexible model. Proposition 1. The multi-variate zero-mean GGD can be decomposed as

G(z; 0 P , Σ, ν) = P j=1 G((U t z) j ; 0, λ j , ν j ) (11)
where, for x = U t z, the distribution of each of its components is given as

G(x; 0, λ, ν) = κ 2λ ν exp - |x| λ ν ν where κ = ν Γ(1/ν) and λ ν = λ Γ(1/ν) Γ(3/ν) ,
where x → G(x; 0, λ, ν) is a real, even, unimodal, bounded and continuous probability density function. It is also differentiable everywhere except for x = 0 when ν 1.

The proof follows directly by injecting the eigen decomposition of Σ in ( 9) and basic properties of x → |x| ν . Proposition 1 shows that, for each of the K clusters, eq. ( 9), indeed, models a prior that is separable in a coordinate system obtained by applying the whitening transform U t . Not only is the prior separable for each coordinate j, but the shape (ν j ) and scale (λ j ) of the distribution may vary.

Before detailing the usage of GGMM priors in EPLL framework, we digress briefly to explain the procedure we used for training such a mixture of generalized Gaussian distributions where different scale and shape parameters are learned for each feature dimension.

Learning GGMMs.

Parameter estimation is carried out using a modified version of the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. EM is an iterative algorithm that performs at each iteration two steps, namely Expectation step (E-Step) and Maximization stop (M-Step), and is known to monotonically increase the model likelihood and converge to a local optimum. For applying EM to learn a GGMM, we leverage standard strategies used for parameter estimation for GGD and/or GGMM that are reported in previous works [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Birney | On the modeling of dct and subband image data for compression[END_REF][START_REF] Sharifi | Estimation of shape parameter for generalized gaussian distributions in subband decompositions of video[END_REF][START_REF] Aiazzi | Estimation based on entropy matching for generalized gaussian pdf modeling[END_REF][START_REF] Dominguez-Molina | A practical procedure to estimate the shape parameter in the generalized gaussian distribution[END_REF][START_REF] Kokkinakis | Exponent parameter estimation for generalized gaussian probability density functions with application to speech modeling[END_REF][START_REF] Boubchir | Multivariate statistical modeling of images with the curvelet transform[END_REF][START_REF] Mohamed | Generalized gaussian mixture model[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF][START_REF] Krupiński | Approximated fast estimator for the shape parameter of generalized gaussian distribution for a small sample size[END_REF]. Our M-Step update for the shape parameter ν is inspired from Mallat's strategy [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] using statistics of the first absolute and second moments of GGDs. Since this strategy uses the method of moments for ν instead of maximum likelihood estimation, we refer to this algorithm as modified EM and the M-Step as Moment step. We also noticed that shape parameters ν < .3 lead to numerical issues and ν > 2 leads to local minima with several degenerate components. For this reason, at each step, we impose the constraint that the learned shape parameters satisfy ν ∈ [. [START_REF] Aiazzi | Estimation based on entropy matching for generalized gaussian pdf modeling[END_REF][START_REF] Aharon | k-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. This observation is consistent with earlier works that have attempted to learn GGMM shape parameters from data [START_REF] Roenko | Estimation of parameters for generalized gaussian distribution[END_REF]. Given n training clean patches of size P and an initialization for the K parameters w k > 0, Σ k ∈ R P ×P and ν k ∈ R P , for k = 1, . . . , K, our modified EM algorithm iteratively alternates between the following two steps:

• Expectation step (E-Step)
-For all components k = 1, . . . , K and training samples i = 1, . . . , n, compute:

ξ k,i ← w k G(z i ; 0 P , Σ k , ν k ) K l=1 w l G(z i ; 0 P , Σ l , ν l ) . • Moment step (M-Step)
-For all components k = 1, . . . , K, update:

w k ← n i=1 ξ k,i K l=1 n i=1 ξ l,i and Σ k ← n i=1 ξ k,i z i z t i n i=1 ξ k,i .
-Perform eigen decomposition of Σ k :

Σ k = U k Λ k U t k where Λ k = diag(λ k,1 , λ k,2 , . . . , λ k,P ) 2 .
-For all components k = 1, . . . , K and dimensions j = 1, . . . , P , compute:

χ k,j ← n i=1 ξ k,i |(U t k z i ) j | n i=1 ξ k,i and (ν k ) j ← Π [.3,2] F -1 χ 2 k,j λ 2 k,j
.

where

Π [a,b] [x] = min(max(x, a), b) and F (x) = Γ(2/x) 2 Γ(3/x)Γ(1/x
) is a monotonic invertible function that was introduced in [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] (we used a lookup table to perform its inversion as done in [START_REF] Birney | On the modeling of dct and subband image data for compression[END_REF][START_REF] Sharifi | Estimation of shape parameter for generalized gaussian distributions in subband decompositions of video[END_REF]). Note that χ 2 k,j and λ 2 k,j corresponds to the first absolute and second moments for component k and dimension j, respectively.

For consistency purposes, we keep the training data and the number of mixture components in the models the same as that used in the original EPLL algorithm [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. Specifically, we train our models on n = 2 million clean patches randomly extracted from Berkeley Segmentation Dataset (BSDS) [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. We learn K = 200 zero-mean generalized Gaussian mixture components from patches of size P = 8 × 8. We opted for a warm-start training by initializing our GGMM model with the GMM model from [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] and with initial values of shape parameters as 2. We run our modified EM algorithm for 100 iterations. As observed in Figure 1, the obtained GGMM models the underlying distributions of a cluster of clean patches much better than a GMM. In addition, we will see in Section 6 that our GGMM estimation did not lead to overfitting as it is also a better fit than a GMM for unseen clean patches.

Patch denoising with GGMM priors.

We now explain why solving step [START_REF] Birney | On the modeling of dct and subband image data for compression[END_REF] in EPLL is non-trivial when using a GGMM patch prior. In this case, for a noisy patch z with variance σ 2 , equation ( 5) becomes ẑ ← argmin

z∈R P 1 2σ 2 ||z -z|| 2 -log K k=1 w k G(z; 0 P , Σ k , ν k ) . ( 12
)
As for GMMs, due to the multi-modality of the GGMM prior, this optimization problem is highly non-convex. To circumvent this issue, we follow the strategy used by EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] in the specific case of Gaussian mixture model prior. The idea is to restrict the sum involved in the logarithm in eq. ( 12) to only one component k .

If we consider the best k to be given (the strategy to select the best k will be discussed next), then eq. ( 12) is approximated by the following simpler problem (13) ẑ ← argmin

z∈R P ||z -z|| 2 2σ 2 -log G(z; 0 P , Σ k , ν k ) = ||z -z|| 2 2σ 2 + ||Σ -1/2 ν k z|| ν ν .
The main advantage of this simplified version is that, by virtue of Proposition 1, the underlying optimization becomes tractable and can be separated into P one-dimensional optimization problems, as:

ẑ = U k x where xj = s(x j ; σ, λ k ,j , ν k ,j ) with x = U t k z (14) and s(x; σ, λ, ν) = s ν σ,λ (x) ∈ argmin t∈R 1 2σ 2 (x -t) 2 + |t| ν λ ν ν where λ ν = λ Γ(1/ν) Γ(3/ν) , (15)
where for all k, ν k,j = (ν k ) j and λ k,j = (λ k ) j . While the problem is not necessarily convex, its solution s ν σ,λ is always uniquely defined almost everywhere (see, Section 5). We call this almost everywhere real function s ν σ,λ : R → R shrinkage function. When ν = 2, it is a linear function that is often referred to as Wiener shrinkage. When ν = 2, as we will discuss in Section 5, it is a non-linear shrinkage function that can be computed in closed form for some cases or with some approximations. Now, we address the question of finding a strategy for choosing a relevant component k to replace the mixture distribution inside the logarithm. The optimal component k can be obtained by maximizing the posterior as

k ∈ argmax 1 k K p(k | z) = argmax 1 k K w k p(z | k) = argmin 1 k K -log w k -log p(z | k) (16)
where the weights of the GGMM corresponds to the prior probability w k = p(k). We next use the fact that the patch z (conditioned on k) can be expressed as z = z + n where z and n are two independent random variables from distributions G(0 P , Σ k , ν k ) and N (0 P , σ 2 Id P ) respectively. It follows that the distribution of z is the convolution of these latter two, and then

-log p(z | k) = -log R P G(z -z; 0 P , Σ k , ν k ) • N (z; 0 P , σ 2 Id P ) dz . ( 17
)
We next use Proposition 1 to separate this integration problem into P one-dimensional integration problems. We obtain

-log p(z | k) = P j=1 f ((U t k z) j ; σ, λ k,j , ν k,j ) (18)
where, for x = U t k z, the integration problem of each of its components reads as

f (x; σ, λ, ν) = f ν σ,λ (x) = -log R G(x -t; 0, λ, ν) • N (t; 0, σ 2 ) dt . ( 19
)
We call the real function f ν σ,λ : R → R the discrepancy function which measures the goodness of fit of a GGD to the noisy value x. When ν = 2, this function is quadratic with x. For ν = 2, as we will discuss in Section 4, it is a non-quadratic function, that can be efficiently approximated based on an in-depth analysis of its asymptotic behavior.

Figure 3. Illustration of our extension of EPLL to GGMM priors. The general procedure, illustrated in the top row, is similar to the original EPLL scheme described in Figure 2 but relies on generalized Gaussian distributions instead of Gaussian distributions. The shape of the discrepancy function, illustrated in the second row, depends on the given scale and shape parameters (λ and ν) of the GGD components. In Section 4, we will see that it can be approximated based on six parameters, four of them retrieved from lookup tables (LUTs). Finally, the shrinkage function, illustrated in the bottom row, can be non-linear and depends on the selected GGD component. In Section 5, we will see that it can be approximated by one of five predefined parametric functions depending on the range in which the scale parameter ν lies. The values ν = 1.2 and λ/σ = .25, shown in the bottom row, were chosen for the sake of illustration.

Figure 3 illustrates the details of the patch denoising step under the GGMM-EPLL framework. It shows that the method relies on fast approximations f ν σ,λ and ŝν σ,λ of the discrepancy and shrinkage functions, respectively.

The next two sections are dedicated to the analysis and approximations of the discrepancy function f ν σ,λ and the shrinkage function s ν σ,λ , respectively. 4. Discrepancy function: analysis and approximations. From its definition given in eq. ( 19), the discrepancy function reads for ν > 0, σ > 0 and λ > 0, as

f ν σ,λ (x) = -log 1 √ 2πσ ν 2λ ν Γ(1/ν) -log ∞ -∞ exp - (x -t) 2 2σ 2 exp - |t| λ ν ν dt . ( 20
)
It corresponds to the negative logarithm of the distribution of the sum of a zero-mean generalized Gaussian and a zero-mean Gaussian random variables. When ν = 2, the generalized Gaussian random variable becomes Gaussian, and the resulting distribution is also Gaussian with zero-mean and variance σ 2 + λ 2 , and then

f 2 σ,λ (x) = 1 2 log 2π + log(σ 2 + λ 2 ) + x 2 σ 2 + λ 2 . (21) Remark 1. For ν = 2, a direct consequence of (21) is that -log p(z | k) is
as an affine function of the Mahanalobis distance between z and 0 P for the covariance matrix Σ k + σ 2 Id P :

-log p(z | k) = 1 2 P log 2π + log |Σ k + σ 2 Id P | + zt (Σ k + σ 2 Id P ) -1 z . (22)
When ν = 1, the generalized Gaussian random variable becomes Laplacian, and the distribution resulting from the convolution also has a closed form which leads to the following discrepancy function

f 1 σ,λ (x) = log(2 √ 2λ) - σ 2 λ 2 -log e √ 2x λ erfc x √ 2σ + σ λ + e - √ 2x λ erfc - x √ 2σ + σ λ , (23)
refer to Appendix A for derivation (note that this expression is given in [START_REF] Gazor | A soft voice activity detector based on a laplacian-gaussian model[END_REF]).

To the best of our knowledge, there are no simple expressions for other values of ν. One solution proposed by [START_REF] Soury | New results on the sum of two generalized gaussian random variables[END_REF] is to express this in terms of the bi-variate Fox-H function [START_REF] Mittal | An integral involving generalized function of two variables[END_REF]. This, rather cumbersome expression, is computationally demanding. In practice, this special function requires numerical integration techniques over complex lines [START_REF] Peppas | A new formula for the average bit error probability of dual-hop amplify-and-forward relaying systems over generalized shadowed fading channels[END_REF], and is thus difficult to numerically evaluate it efficiently. Since, in our application, we need to evaluate this function a large number of times, we cannot utilize this solution.

In [START_REF] Soury | New results on the sum of two generalized gaussian random variables[END_REF], the authors have also proposed to approximate this non-trivial distribution by another GGD. For fixed values of σ, λ and ν, they proposed three different numerical techniques to estimate its parameters λ and ν that best approximate either the kurtosis, the tail or the cumulative distribution function. Based on their approach, the discrepancy function f ν σ,λ (x) would thus be a power function of the form |x| ν .

In this paper, we show that f ν σ,λ does, indeed, asymptotically behave as a power function for small and large values of x, but the exponent can be quite different for these two asymptotics. We believe that these different behaviors are important to be preserved in our application context. For this reason, f ν σ,λ cannot be modeled as a power function through a GGD distribution. Instead, we provide an alternative solution that is able to capture the correct behavior for both of these asymptotics, and that also permits fast computation.

4.1. Theoretical analysis. In this section, we perform a thorough theoretical analysis of the discrepancy function, in order to approximate it accurately. Let us first introduce some basic properties regarding the discrepancy function.

Proposition 2. Let ν > 0, σ > 0, λ > 0 and f ν σ,λ as defined in eq. [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. The following relations hold true

f ν σ,λ (x) = log σ + f ν 1,λ/σ (x/σ) , (reduction) f ν σ,λ (x) = f ν σ,λ (-x) , (even) |x| |y| ⇔ f ν σ,λ (|x|) f ν σ,λ (|y|) , (unimodality) min x∈R f ν σ,λ (x) = f ν σ,λ (0) > -∞ . (lower bound at 0)
The proofs can be found in Appendix B. Based on Proposition 2, we can now express the discrepancy function f ν σ,λ (x) : R → R in terms of a constant γ ν λ and another function ϕ ν λ : R * + → R, both of which can be parameterized by only two parameters λ > 0 and ν > 0, as

f ν σ,λ (x) = log σ + γ ν λ/σ + e ϕ ν λ/σ (|x/σ|) if x = 0 , 0 otherwise , (24)
where

ϕ ν λ (x) = log f ν 1,λ (x) -γ ν λ and γ ν λ = f ν 1,λ (0) . ( 25
)
We call ϕ ν λ the log-discrepancy function. At this point, let us consider an instructive toy example for the case when ν = 2. In this case, from eq. ( 21), we can deduce that the log-discrepancy function is a log-linear function (i.e., a linear function of log x)

ϕ 2 λ (x) = α log x + β , (26)
and

γ 2 λ = 1 2 log 2π + log(1 + λ 2 ) , (27)
where α = 2 and β = -log 2 -log(1 + λ 2 ) . [START_REF] Huber | Robust estimation of a location parameter[END_REF] Here, the slope α = 2 reveals the quadratic behavior of the discrepancy function. Figure 4 gives an illustration of the resulting convolution (a Gaussian distribution), the discrepancy function (a quadratic function) and the log-discrepancy (a linear function with slope 2). Note that quadratic metrics are well-known to be non-robust to outliers, which is in complete agreement with the fact that Gaussian priors have thin tails.

Another example is the case of ν = 1. From eq. ( 23), the log-discrepancy is given by and

ϕ 1 λ (x) = log log 2 erfc 1 λ -log e √ 2x λ erfc x √ 2 + 1 λ + e - √ 2x λ erfc -x √ 2 + 1 λ , (29)
γ 1 λ = 1 2 log 2 + log λ - 1 λ 2 -log erfc 1 λ . ( 30
)
Unlike for ν = 2, this function is not log-linear and thus f 1 σ,λ is not a power function. Nevertheless, as shown by the next two theorems, it is also asymptotically log-linear for small and large values of x.

Theorem 1. The function ϕ 1

λ is asymptotically log-linear in the vicinity of 0

ϕ 1 λ (x) ∼ 0 α 1 log x + β 1 , (31)
where α 1 = 2 and

β 1 = -log λ + log 1 √ π exp -1 λ 2 erfc 1 λ - 1 λ . ( 32
)
The proof can be found in Appendix C.

Theorem 2. The function ϕ 1

λ is asymptotically log-linear in the vicinity of +∞

ϕ 1 λ (x) ∼ ∞ α 2 log x + β 2 , ( 33
)
where α 2 = 1 and

β 2 = 1 2 log 2 -log λ . (34)
The proof can be found in Appendix D.

Theorem 1 and Theorem 2 show that ϕ 1 λ has two different asymptotics that can be approximated by a log-linear function. Interestingly, the exponent α 1 = 2 in the vicinity of 0 shows that the Gaussian distribution involved in the convolution prevails over the Laplacian distribution and thus, the behavior of f 1 σ,λ is quadratic. Similarly, the exponent α 2 = 1 in the vicinity of +∞ shows that the Laplacian distribution involved in the convolution prevails over the Gaussian distribution and the behavior of f 1 σ,λ is then linear. These results are supported by Figure 5 which illustrates the resulting convolution, the discrepancy function (eq. ( 23)) and the log-discrepancy function (eq. ( 29)). Furthermore, the discrepancy function f 1 σ,λ shares a similar behavior with the well-known Huber loss function [START_REF] Huber | Robust estimation of a location parameter[END_REF] (also called smoothed 1 ), known to be more robust to outliers. This is again in complete agreement with the fact that Laplacian priors have heavier tails.

In the case 2 3 < ν < 2, even though ϕ ν λ has no simple closed form expression, the similar conclusions can be made as a result of the next two theorems. Theorem 3. Let ν > 0. The function ϕ ν λ is asymptotically log-linear in the vicinity of 0

ϕ ν λ (x) ∼ 0 α 1 log x + β 1 ,
where α 1 = 2 and β 1 = -log 2 + log

    1 - ∞ -∞ t 2 e -t 2 2 exp - |t| λ ν ν dt ∞ -∞ e -t 2 2 exp - |t| λ ν ν dt     .
The proof can be found in Appendix E.

Theorem 4. Let 2 3 < ν < 2, then ϕ ν λ is asymptotically log-linear in the vicinity of +∞

ϕ ν λ (x) ∼ ∞ α 2 log x + β 2 ,
where α 2 = ν and

β 2 = -ν log λ - ν 2 log Γ(1/ν) Γ(3/ν) .
The proof relies on a result of Berman (1992) [START_REF] Berman | The tail of the convolution of densities and its application to a model of HIV-latency time[END_REF] and is detailed in Appendix F.

Remark 2. For ν > 2, an asymptotic log-linear behavior with α 2 = 2 and β 2 = -log 2 can be obtained using exactly the same sketch of proof as the one of Theorem 4. Remark 5. For 0 < ν Again, the exponent α 1 = 2 in the vicinity of 0 shows that the Gaussian distribution involved in the convolution prevails over the generalized Gaussian distribution and the behavior of f ν σ,λ is then quadratic. Similarly, the exponent α 2 = ν in the vicinity of +∞ shows that the generalized Gaussian distribution involved in the convolution prevails over the Gaussian distribution and the behavior of f ν σ,λ is then a power function of the form x ν . These results are supported by Figure 6 that illustrates the resulting convolution, the discrepancy function and the log-discrepancy function for ν = 1.5 and ν = .3. Moreover, the discrepancy function f ν σ,λ with ν 1 shares a similar behavior with well-known robust M-estimator loss functions [START_REF] Huber | Robust statistics[END_REF]. In particular, the asymptotic case for ν → 0 resembles the Tukey's bisquare loss, known to be insensitive to outliers. This is again in complete agreement with GGD priors having larger tails as ν goes to 0. Figure 7 shows the evolution of the log-discrepancy function for various values of ν in the context of three different signal-to-noise ratios λ/σ (SNR). One can observe that as the SNR decreases (resp., increases), the left (resp., right) asymptotic behavior starts dominating over the right (resp., left) asymptotes. In other words, for ν < 2, the intersection of the two asymptotes goes to +∞ (resp., -∞). Last but not least, for 0 < ν 2, the log-discrepancy function ϕ ν λ is always concave and since α 2 α 1 it is thus upper-bounded by its left and right asymptotes.

From Theorem 3, Theorem 4 and Remark 5, we can now build two asymptotic log-linear approximations for ϕ ν λ , with 0 < ν 2, and subsequently an asymptotic power approximation for f ν σ,λ by using the relation [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF]. Next, we explain the approximation process of the inbetween behavior, as well as its efficient evaluation.

Numerical approximation.

We now describe the proposed approximation of the discrepancy function f ν 1,λ through an approximation φν λ of the log-discrepancy function as

f ν 1,λ (x) = γ ν λ + exp φν λ (x) where γ ν λ = f ν 1,λ (0) . (35)
Based on our previous theoretical analysis, a solution preserving the asymptotic, increasing and concave behaviors of ϕ ν 1,λ can be defined by making use of the following approximations

φν λ (x) = α 1 log |x| + β 1 -rec(α 1 log |x| + β 1 -α 2 log |x| -β 2 ) , ( 36
)
where rec is a so-called rectifier function that is positive, increasing, convex and satisfies lim x→-∞ rec(x) = 0 and rec(x) ∼ x→∞ x . [START_REF] Michaeli | Blind deblurring using internal patch recurrence[END_REF] In this paper, we consider the two following rectifying functions relu(x) = max(0, x) and softplus(x . Illustrations of our approximations of ϕ ν λ and the corresponding underlying posterior distribution N (0, ν, λ) * G(0, 1) (where ν = .8 and λ = 4). The blue curves have been obtained by evaluating the convolution using numerical integration techniques for all x. The dashed curves are obtained using the proposed relu-and softplus-based approximations that have closed-form expressions.

) = h log 1 + exp x h , h > 0 , (38)
as coined respectively in [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] and [START_REF] Dugas | Incorporating second-order functional knowledge for better option pricing[END_REF]. Using the function relu (Rectified linear unit) leads to an approximation φν λ that is exactly equal to the asymptotes of ϕ ν λ with a singularity at their crossing point. In this paper, we will instead use the function softplus as it allows the approximation of ϕ ν λ to converge smoothly to the asymptotes without singularity. Its behavior is controlled by the parameter h > 0. The smaller the value of h is, the faster the convergence speed to the asymptotes.

The parameter h should be chosen such that the approximation error between φν λ (x) and ϕ ν λ (x) is as small as possible. This can be done numerically by first evaluating ϕ ν λ (x) with integration techniques for a large range of values x, and then selecting the parameter h by least square. Of course, the optimal value for h depends on the parameter λ and ν.

Figure 8 gives an illustration of our approximations of the log-discrepancy and the corresponding distribution obtained with relu and softplus. On this figure the underlying functions have been obtained by numerical integration for a large range of value of x. One can observe that using softplus provides a better approximation than relu.

Our approximation for f ν 1,λ (x) is parameterized by six scalar values: γ λ ν , α 1 , β 1 , α 2 , β 2 and h that depend only on the original parameters λ and ν. From our previous analysis, we have that α 1 = 2 and α 2 = ν. The other parameters are non-linear functions of λ and ν. The parameters γ λ ν , β 1 and β 2 require either performing numerical integration or evaluating the special function Γ. As discussed, the parameter h requires numerical integration for various x and then optimization. For these reasons, these values cannot be computed during runtime. Instead, we pre-compute these four parameters offline for 10, 000 different combinations of λ and ν values in the intervals [10 -3 , 10 3 ] and [0.3, 2], respectively (the choice for this range was motivated in Subsection 3.1). The resulting values are then stored in four corresponding lookup tables. During runtime, these parameters are retrieved online by bi-linear extrapolation and interpolation. The four lookup tables are given in Figure 9. We will see in Section 6 that using the approximation f ν 1,λ results in substantial acceleration without significant loss of performance as compared to computing f ν 1,λ directly by numerical integration during runtime.

Shrinkage functions: analysis and approximations.

Recall that from its definition given in eq. ( 15), the shrinkage function is defined for ν > 0, σ > 0 and λ > 0, as 5.1. Theoretical analysis. Except for some particular values of ν (see, Subsection 5.2), Problem [START_REF] Mohamed | Generalized gaussian mixture model[END_REF] does not have explicit solutions. Nevertheless, as shown in [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF], Problem (39) admits two (not necessarily distinct) solutions. One of them is implicitly characterized as

s ν σ,λ (x) ∈ argmin t∈R (x -t) 2 2σ 2 + λ -ν ν |t| ν . (39)
s ν σ,λ (x) = 0 if 0 < ν 1 and |x| τ ν λ , t otherwise , ( 40
)
where t = x -sign(t )νσ 2 λ -ν ν |t | ν-1 ,
and

τ ν λ = (2 -ν)(2 -2ν) -1-ν 2-ν (σ 2 λ -ν ν ) 1 2-ν if ν < 1 , σ 2 λ -1 otherwise (ν = 1) .
The other one is obtained by changing |x| τ ν λ to |x| < τ ν λ in [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF], and so they coincide almost every (x, λ, σ, ν). As discussed in [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF], for ν > 1, s ν σ,λ (x) is differentiable, and for ν 1, the shrinkage exhibits a threshold τ ν λ that produces sparse solutions. Proposition 3 summarizes a few important properties. Proposition 3. Let ν > 0, σ > 0, λ > 0 and s ν σ,λ as defined in eq. [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian density and kullback-leibler distance[END_REF]. The following relations hold true

s ν σ,λ (x) = σs ν 1, λ σ x σ , (reduction)
s ν σ,λ (x) = -s ν σ,λ (-x) , (odd) s ν σ,λ (x) ∈ [0, x] if x 0 [x, 0] otherwise , (shrinkage) x 1 x 2 ⇔ s ν σ,λ (x 1 ) s ν σ,λ (x 2 ) , (increasing with x) λ 1 λ 2 ⇔ s ν σ,λ 1 (x) s ν σ,λ 2 (x) , (increasing with λ) lim λ σ →0 s ν σ,λ (x) = 0 , (kill low SNR)
ν Shrinkage s ν σ,λ (x) Remark < 1 x -γx ν-1 + O(x 2(ν-1) ) if |x| τ ν λ 0 otherwise ≈ Hard-thresholding [40] 1 sign(x) max |x| - √ 2σ 2 λ , 0 Soft-thresholding [17] 4/3 x + γ 3 ζ -x 2 - 3 ζ + x 2 [10]
3/2 sign(x)

γ 2 + 4|x| -γ 2 4 [10] 2 λ 2 λ 2 + σ 2 • x Wiener (LMMSE) with γ = νσ 2 λ -ν ν and ζ = x 2 + 4 γ 3 3 . lim λ σ →+∞ s ν σ,λ (x) = x . (keep high SNR)
The proofs can be found in Appendix G. These properties show that s ν σ,λ is indeed a shrinkage function (non-expansive). It shrinks the input coefficient x according to the model ν and the modeled signal to noise ratio λ σ (SNR). When x is small in comparison to the SNR, it is likely that its noise component dominates the underlying signal, and is, therefore, shrunk towards 0. Similarly, when x is large, it will likely be preserved. This is even more likely when ν is small, since in this case large coefficients are favored by the prior. Illustrations of shrinkage functions for various SNR and ν are given in Figure 10.

Numerical approximations. The shrinkage function s ν

σ,λ , implicitly defined in (40) does not have a closed form expression in general. Nevertheless, for fixed values of x, σ, λ, ν, s ν σ,λ (x) can be estimated using solvers such as Newton descent or Halleys root-finding method. These approaches converge quite fast and, in practice, reach a satisfying solution within ten iterations. However, since in our application of interest we need to evaluate this function a large number times, we will follow a different path in order to reduce computation time (even though we have implemented this strategy).

As discussed earlier, s ν σ,λ is known in closed form for some values of ν, more precisely: ν = {1, 4/3, 3/2, 2} (as well as ν = 3 but this is out of the scope of this study), see for instance [START_REF] Chaux | A variational formulation for framebased inverse problems[END_REF]. When ν = 2, we retrieve the linear minimum mean square estimator (known in signal Figure 11. Set of 100 patches, sorted by the norm of their gradient, and generated to be independently distributed according to (from left to right) a GMM, GGMM, LMM and HLMM. For ease of visualization, only the top eigendirections corresponding to 80% of the variance have been chosen. Near-constant patches with variance smaller than 2 P ||Σ k || 2 F have also been discarded.

processing as Wiener filtering) and related to Tikhonov regularization and ridge regression. This shrinkage is linear and the slope of the shrinkage goes from 0 to 1 as the SNR increases (see Figure 10). When ν = 1, the shrinkage is the well-known soft-thresholding [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], corresponding to the maximum a posteriori estimator under a Laplacian prior. When ν < 1, the authors of [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF] have shown that (i) the shrinkage admits a threshold with a closed-form expression (given in eq. ( 40)), and (ii) the shrinkage is approximately equal to hard-thresholding with an error term that vanishes when |x| → ∞. All these expressions are summarized in Table 1.

In order to keep our algorithm as fast as possible, we propose to use the approximation of the shrinkage given for ν < 1 in [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF]. Otherwise, we pick one of the four shrinkage functions corresponding to ν = {1, 4/3, 3/2, 2} by nearest neighbor on the actual value of ν ∈ [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Aharon | k-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. Though this approximation may seem coarse compared to the one based on iterative solvers, we did not observe any significant loss of quality in our numerical experiments (see Section 6). Nonetheless, this alternative leads to 6 times speed-up while evaluating shrinkage.

6. Experimental evaluation. In this section we explain the methodology used to evaluate the GGMM model, and present numerical experiments to compare the performance of the proposed GGMM model over existing GMM-based image denoising algorithms. To demonstrate the advantage of allowing for a flexible GGMM model, we also present results using GGMM models with fixed shape parameters, ν = 1 (Laplacian mixture model) and ν = 0.5 (Hyper-Laplacian mixture model). For learning Laplacian mixture model (LMM) and hyper-Laplacian mixture model (HLMM), we use the same procedure as described in Subsection 3.1 but force all shape parameters to be equal to 1 or 0.5, respectively.

Model validation. As discussed in Section 3, Figure 1 illustrates the validity of our model choices with histograms of different dimensions of a single patch cluster. It clearly shows the importance of allowing the shape and scale parameter to vary across dimensions for capturing underlying patch distributions. Since GGMM (and obviously, GMM) falls into the class of generative models, one can also assess the expressivity of a model by analyzing the variability of generated patches and its ability to generate relevant image features (edges, texture elements etc.). This can be tested by selecting a component k of the GGMM (or GMM) with probability w k and sampling patches from the GGD (or GD) as described in [START_REF] Nardon | Simulation techniques for generalized gaussian densities[END_REF]. Figure 11 presents a collage of 100 patches independently generated by this procedure using GMM, GGMM, LMM and HLMM. As observed, patches generated from GGMM show greater balance between strong/faint edges, constant patches and subtle textures than the models that use constant shape parameters such as GGM, LMM and HLMM.

The superiority of our GGMM model over GMM, LMM or HLMM models can also be illustrated by comparing the log-likelihood (LL) achieved by these models over a set of clean patches from natural images. Note that, to maintain objectivity, the models have to be tested on data that is different than the dataset used during training. To this end, we compute the LL of the four above-mentioned models on all non-overlapping patches of 40 randomly selected images extracted from BSDS testing set [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF], which is a different set than the training images used in the EM algorithm (parameter estimation/model learning). One can observe that not only GGMM is a better fit than GMM, LMM and HLMM on average for the 40 images, but it is also a better fit on each single image. Given that GGMM have a larger degree of freedom than GMM, this study proves that our learning procedure did not fall prey to over-fitting, and that the extra flexibility provided by GGMM was used to capture relevant and accurate image patterns.

Denoising evaluation. Following the verification of the model, we provide a thorough evaluation of our GGMM prior in denoising task by comparing its performance against EPLL that uses a GMM prior [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] and with our LMM and HLMM models explained above. For the ease of comparison, we utilize the pipeline and settings that was prescribed for the original EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] algorithm (see Section 2.1). To reduce the computation time of all EPLL-based algorithms, we utilize the random patch overlap procedure introduced by [START_REF] Parameswaran | Accelerating gmm-based patch priors for image restoration: Three ingredients for a 100× speed-up[END_REF]. That is, instead of extracting all patches at each iteration, a randomly selected but different subset of overlapping patches consisting of only 3% of all possible patches is processed in each iteration. For the sake of reproducibility of our results, we have made our MATLAB/MEX-C implementation available online at https://bitbucket.org/cdeledalle/ggmm-epll. The evaluation is carried out on classical images such as Barbara, Cameraman, Hill, House, Lena and Mandrill, and on 60 images taken from BSDS testing set [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] (the original BSDS test data contains 100 images, the other 40 were used for model validation experiments). All image have been corrupted independently by ten independent realizations of additive white Gaussian noise with standard deviation σ = 5, 10, 20, 40 and 60 (with pixel values between [0, 255]). The EPLL algorithm using mixture of Gaussian, generalized Gaussian priors are indicated as GMM and GGMM in Table 2. Results obtained with BM3D algorithm [START_REF] Dabov | Image denoising by sparse 3-d transformdomain collaborative filtering[END_REF] are also included for reference purposes. To stay with the focus of this paper, i.e., on the effect of image priors on EPLL-based algorithms, BM3D will be excluded from our performance comparison discussions. The denoising performance of the algorithms are measured in terms of Peak Signal to Noise Ratio (PSNR) and Structural SIMilarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. As can be observed in Table 2, in general, GGMM prior provides better PSNR performance than the GMM prior. In terms of SSIM values (shown in the bottom part of Table 2), GGMM prior is comparable to GMM. In order to demonstrate the effect of fixed ν values compared to the more flexible GGMM prior, we compare the results of GGMM against GMM (ν = 2), Laplacian Mixture Model (GGMM with ν = 1) and hyper-Laplacian mixture model (GGMM with ν = 0.5) priors in the same scenarios for σ = 5, 20 and 60. These results are shown in Table 3. GGMM prior provides better PSNR performance on average than the fixed-shape priors. The differences in denoising performance can also be verified visually in Figure 13, Figure 14 and Figure 15. The denoised images obtained using GGMM prior show much fewer artifacts as compared to GMM-EPLL results, in particular in homogeneous regions. On the other hand, GGMM prior is also able to better preserve textures than LMM and HLMM.

Prior fitness for image denoising. In this work, we have considered non-blind image denoising. That is, the noise standard deviation is assumed to be known. In this setting, if the restoration model is accurate, one should ideally achieve optimal restoration performance when using the true degradation. To verify this, we conducted a denoising task with image corrupted with noise with standard deviation σ = 20. We used GMM, LMM and GGMM priors in the restoration framework with assumed σ values ranging from 15 to 30. Figure 16 shows the evolution of average restoration performance over 40 images from BSDS testing set (kept aside for validation, as mentioned above) with varying noise variances. GGMM prior attains its best performance when the noise variance used in the restoration model matches with the ground truth σ = 20. In contrast to GGMM, GMM (resp., LMM) reaches its best performance at a larger (resp., lower) value of σ than the correct noise used during degradation. This is because GMM tends to under-smooth clean patches (resp., over-smooth) so that a larger (resp., lower) value of σ is required to compensate the mismatch between the assumed prior and the actual distribution in the restoration model. This indicates that GGMM is a better option to model distribution of image patches than GMM or LMM.

Influence of our approximations. All previous experiments using GGMM patch priors were conducted based on the two proposed approximations introduced in Section 4 and Section 5. In Figure 17 and Table 4, we provide a quantitative illustration of the speed-ups provided by these approximations and their effect on denoising performance. Timings were carried out with Matlab 2018a on an Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz (neither multi-core paralellization nor GPUs acceleration were used). Figure 17a shows the result obtained by calculating original discrepancy function via numerical integration and the shrinkage function via Halley's root-finding method. This makes the denoising process extremely slow and takes 10 hours and 29 minutes for denoising an image of size 128×128 pixels. The approximated discrepancy function provides 4 orders of magnitude speed-up with no perceivable drop in performance (Figure 17b). In addition, incorporating the shrinkage approximation provides further acceleration that allows the denoising to complete in less than 2 seconds with a very minor drop in PSNR/SSIM. As indicated in the detailed profiles on Table 4, the shrinkage

Table 4

Runtime profiles (averaged over ten runs) of GGMM-EPLL corresponding to the denoising experiment shown in Figure 17. These profiles are obtained respectively by evaluating the classification and shrinkage problem with numerical solvers (numerical integration and Halley's root-finding method), approximating the classification problem only, and approximating both problems. Profiles are split into discrepancy, shrinkage and patch extration/reprojection. Speed-up with respect to the no-approximation (first column) are indicated in parenthesis and major accelerations in green. Percentage of time taken for each step with respect to the overall execution time (first row), are indicated below each time reading and bottlenecks are indicated in red. approximation provides an acceleration of six-fold to the shrinkage calculation step itself and leads to an overall speed-up of 1.5 due to the larger bottleneck caused by discrepancy function calculation. The approximately 23,000× speed-up realized without any perceivable drop in denoising performance underscores the efficacy of our proposed approximations.

No approximations

7. Conclusions and Discussion. In this work, we suggest using a mixture of generalized Gaussians for modeling the patch distribution of clean images. We provide a detailed study of the challenges that one encounters when using a highly flexible GGMM prior for image restoration in place of a more common GMM prior. We identify the two main bottlenecks in the restoration procedure when using EPLL and GGMM -namely, the patch classification step and the shrinkage step. One of the main contributions of this paper, is the thorough theoretical analysis of the classification problem allowing us to introduce an asymptotically accurate approximation that is computationally efficient. In order to tackle the shrinkage step, we collate and extend the existing solutions under GGMM prior.

Our numerical experiments indicate that our flexible GGMM patch prior is a better fit for modeling natural images than GMM and other mixture distributions with constant shape parameters such as LMM or HLMM. In image denoising tasks, we have shown that using GGMM priors, often, outperforms GMM when used in the EPLL framework.

Nevertheless, we believe the performance of GGMM prior in these scenarios falls short of its expected potential. Given that GGMM is persistently a better prior than GMM (in terms of log-likelihood), one would expect the GGMM-EPLL to outperform GMM-EPLL consistently. We postulate that this under-performance is caused by the EPLL strategy that we use for optimization. That is, even though the GGMM prior may be improving the quality of the global solution, the half quadratic splitting strategy used in EPLL is not guaranteed to return a better solution due to the non-convexity of the underlying problem. For this reason, we will focus our future work in designing specific optimization strategies for GGMM-EPLL leveraging the better expressivity of the proposed prior model for denoising and other general restoration applications.

Another direction of future research will focus on extending this work to employ GGMMs/GGDs in other model-based signal processing tasks. Of these tasks, estimating the parameters of GGMM directly on noisy observations is a problem of particular interest, that could benefit from our approximations. Learning GMM priors on noisy patches has been shown to be useful in patch-based image restoration when clean patches are not available a priori, or to further adapt the model to the specificities of a given noisy image [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF][START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising (HDMI)[END_REF].

Another open problem is to analyze the asymptotic behavior of the minimum mean square estimator (MMSE) shrinkage with GGD prior, as an alternative to MAP shrinkage. This could be useful to design accurate approximations for other general inference frameworks. Last but not least, characterizing the exact asymptotic behaviors of the convolution of two arbitrary GGDs, as investigated in [START_REF] Soury | New results on the sum of two generalized gaussian random variables[END_REF], is still an open question. To the best of our knowledge, our study is the first attempt towards this goal but in ours one of the GGD is always Gaussian (noise). Extending our study to the general GGD case (or even specific cases such as Laplacian) is a challenging problem that is of major interest in signal processing tasks where noise is not Gaussian but instead follows another GGD.

Appendix A. Proof of equation [START_REF] Gazor | Speech probability distribution[END_REF].

Proof. For ν = 1, using Γ(1) = 1 and Γ(3) = 2, we obtain

f 1 σ,λ (x) = log(2 √ πσλ) -log ∞ -∞ e -t 2 2σ 2 - √ 2|x-t| λ dt , (41) = log(2 √ πσλ) -log e - √ 2x λ x -∞ e -t 2 2σ 2 + √ 2t λ dt + e √ 2x λ ∞ x e -t 2 2σ 2 - √ 2t λ dt . (42) Since erf (t) = 2e -t 2 √ π , it follows that for a > 0 and b > 0 ∂ ∂t - √ πa 2 e a 4b 2 erf - t √ a - √ a 2b = e -t 2 a -t b . (43)
Therefore we have with a = 2σ 2 and b

= λ/ √ 2 x -∞ e -t 2 2σ 2 + √ 2t λ dt = - √ πσe σ 2 λ 2 √ 2 erf - t √ 2σ + σ λ x -∞ = √ πσe σ 2 λ 2 √ 2 erfc - x √ 2σ + σ λ , (44)
since lim t→∞ erf(t) = 1 and erfc(t) = 1 -erf(t). Similarly, we get

∞ x e -t 2 2σ 2 - √ 2t λ dt = - √ πσe σ 2 λ 2 √ 2 erfc x √ 2σ + σ λ . (45)
Plugging these two last expressions in [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] and rearranging the terms conclude the proof.

Appendix B. Proof of Proposition 2.

Proof. Starting from the definition of f ν σ,λ and using the change of variable t → σt, eq. (reduction) follows as [START_REF] Parameswaran | Accelerating gmm-based patch priors for image restoration: Three ingredients for a 100× speed-up[END_REF]

f ν σ,λ (x) = -log R σ κ 2λ ν exp - σ|t| λ ν ν 1 √ 2πσ exp - (x -σt) 2 2σ 2 dt , = log σ-log R κ 2λ ν /σ exp - |t| λ ν /σ ν 1 √ 2π exp - (x/σ -t) 2 2 dt = log σ+f ν 1,λ/σ x σ .
Properties (even) and (unimodality) hold since the convolution of two real even unimodal distributions is even unimodal [START_REF] Wintner | Asymptotic distributions and infinite convolutions. edwards brothers[END_REF][START_REF] Purkayastha | Simple proofs of two results on convolutions of unimodal distributions[END_REF]. Property (lower bound at 0) follows from (even), (unimodality) and the fact that the convolution of continuous and bounded real functions are continuous and bounded.

Appendix C. Proof of Theorem 1.

Lemma C.1. Let a > 0 and b > 0. For x in the vicinity of 0, we have

1 2abx log erfc(ax + b) + e -4abx erfc(-ax + b) 2 erfc(b) = -1 + ab - ae -b 2 erfc(b) √ π x + o(x) . (47)
Proof. Since erfc (x) = -2e -x 2 √ π and erfc (x) = 2xe -x 2 √ π , using second order Taylor's expansion for x in the vicinity of 0, it follows that

erfc(ax + b) ∼ 0 erfc(b) - 2ae -b 2 √ π x + 2a 2 be -b 2 √ π x 2 , ( 48
) and erfc(-ax + b) ∼ 0 erfc(b) + 2ae -b 2 √ π x + 2a 2 be -b 2 √ π x 2 . ( 49
)
We next make the following deductions

e -4abx erfc(-ax + b)) ∼ 0 e -4abx erfc(b) + e -4abx 2ae -b 2 √ π x + e -4abx 2a 2 be -b 2 √ π x 2 , (50) erfc(ax+b)+e -4abx erfc(-ax+b)) ∼ 0 (1+e -4abx ) erfc(b)+ 2a 2 be -b 2 √ π x 2 -(1-e -4abx ) 2ae -b 2 √ π x , erfc(ax+b)+e -4abx erfc(-ax+b)) 2 erfc b A(x) ∼ 0 1+e -4abx 2 1+ 2a 2 be -b 2 erfc(b) √ π x 2 -(1-e -4abx ) ae -b 2 erfc(b) √ π x .
The left-hand side A(x) of this last equation is then, in the vicinity of x = 0, equals to

A(x) = 1 + e -4abx 2 1 + 2a 2 be -b 2 erfc(b) √ π x 2 -(1 -e -4abx ) ae -b 2 erfc(b) √ π x + o(x 2 ) . (51)
We next use second-order Taylor's expansion for e -4abx , leading to

A(x) = (1 -2abx + 4a 2 b 2 x 2 + o(x 2 )) 1 + 2a 2 be -b 2 erfc(b) √ π x 2 (52) -(4abx -8a 2 b 2 x 2 + o(x 2 )) ae -b 2 erfc(b) √ π x + o(x 2 ) , = 1 -2abx + 4a 2 b 2 - 2a 2 be -b 2 erfc(b) √ π x 2 + o(x 2 ) . (53)
By using the second-order Taylor's expansion of log(1 + x), it follows that log

[A(x)] = -2abx + 4a 2 b 2 - 2a 2 be -b 2 erfc(b) √ π x 2 -2a 2 b 2 x 2 + o(x 2 ) . (54)
Dividing both sides by 2abx then concludes the proof,

1 2abx log erfc(ax + b) + e -4abx erfc(-ax + b)) 2 erfc b = -1 + ab - ae -b 2 erfc(b) √ π x + o(x) . (55)
Proof of Theorem 1. We first rewrite ϕ 1 λ (x) as

ϕ 1 λ (x) = log log 2 erfc 1 λ -log e √ 2x λ erfc x √ 2 + 1 λ + e - √ 2x λ erfc - x √ 2 + 1 λ , = log √ 2x λ + log   -1 - λ √ 2x log   erfc x √ 2 + 1 λ + e -2 √ 2x λ erfc -x √ 2 + 1 λ 2 erfc 1 λ     . (56)
Next, using Lemma C.1 with a = 1/ √ 2 and b = 1/λ, it follows that

ϕ 1 λ (x) = log √ 2x λ + log - 1 √ 2λ - e -1 λ 2 √ 2 erfc( 1 λ ) √ π x + o(x) , (57)
= log x 2 λ + log e -1 λ 2 erfc( 1 λ ) √ π - 1 λ + o(1) , (58)
= log x 2 λ + log 1 √ π e -1 λ 2 erfc( 1 λ ) - 1 λ + o(1) , (59)
where the last equation follows from the first-order Taylor expansion of log(a + x).

where we have used the knowledge that f ∼ g implies that log f ∼ log g.

Proof of Theorem 2. By writing ϕ 1 λ as in eq. ( 56) and using Lemma D.1 with a = 1/ √ 2 and b = 1/λ, it follows that

ϕ 1 λ (x) = log √ 2x λ + log 1 + λ √ 2x log erfc 1 λ + o (1) = log √ 2x λ + o (1) . ν 2λ ν Γ(1/ν) -log ∞ -∞ e -x 2 2 e xt e -t 2 2 exp - |t| λ ν ν dt , (67)
= -log 1 √ 2π ν 2λ ν Γ(1/ν) + x 2 2 -log ∞ -∞ ∞ k=0 (xt) k k! e -t 2 2 exp - |t| λ ν ν dt . (68)
For x in the vicinity of 0, we can consider |x| 1, and then

∞ -∞ ∞ k=0 (xt) k k! e -t 2 2 exp - |t| λ ν ν dt ∞ -∞ ∞ k=0 |t| k k! e -t 2 2 exp - |t| λ ν ν dt , (69) ∞ -∞ e -t 2 2 +|t|- |t| λν ν dt < ∞ . (70)
Then, Fubini's theorem applies and we get Using third-order Taylor expansion of log(1 + x) for x in the vicinity of 0, it follows that We first recall in a lemma, a result extracted from Corollary 3.3 in [START_REF] Berman | The tail of the convolution of densities and its application to a model of HIV-latency time[END_REF].

f ν 1,λ (x) = -log 1 √ 2π ν 2λ ν Γ(1/ν) + x 2 2 -log ∞ k=0 ∞ -∞ (xt) k k! e -t 2 2 exp - |t| λ ν ν dt , (71)
= -log 1 √ 2π ν 2λ ν Γ(1/ν) + x 2 2 -log ∞ k=0 x k k! ∞ -∞ t k e -
f ν 1,λ (x) = γ ν λ + x 2 2 -log   1 + x 2 2 ∞ -∞ t 2 e -
Lemma F.1 (Berman). Let p and q be differentiable real probability density functions. Define for x large enough u(x) = p -1 (q(x)) and define v and w as p(x -t)q(t) dt . (83) with q(x) = G(x; 0, λ, ν) and p(x) = N (x; 0, 1). We have We now prove eq. (shrinkage). Let t = s ν σ,λ (x), and since t minimizes the objective, then

λ -ν ν |t| ν (x -t) 2 2σ 2 + λ -ν ν |t| ν (x -x) 2 2σ 2 + λ -ν ν |x| ν = λ -ν ν |x| ν , (96)
which implies that |t| |x|. Let x > 0 and assume t = s ν σ,λ (x) < 0. Since t minimizes the objective, then (x -t) 2 2σ 2 + λ -ν ν |t| ν (x + t) 2 2σ 2 + λ -ν ν |t| ν (97) which implies that -xt xt and leads to a contradiction. Then for x > 0, s ν σ,λ (x) ∈ [0, x], which concludes the proof since s ν σ,λ is odd. We now prove (increasing with x). Let x 1 > x 2 and define t 1 = s ν σ,λ (x 1 ) and t 2 = s ν σ,λ (x 2 ). Since t 1 and t 2 minimize their respective objectives, the following two statements hold

(x 1 -t 1 ) 2 2σ 2 + λ -ν ν |t 1 | ν (x 1 -t 2 ) 2 2σ 2 + λ -ν ν |t 2 | ν , (98)
⇒ t 1 (x 1 -x 2 ) t 2 (x 1 -x 2 ) ⇒ t 1 t 2 (since x 1 > x 2 ) . (102)
We now prove (increasing with λ). Let λ 1 > λ 2 and define t 1 = s ν σ,λ 1 (x) and t 2 = s ν σ,λ 2 (x). Since t 1 and t 2 minimize their respective objectives, the following expressions hold (x -t 1 ) which is in contradiction with the fact that t minimizes the objective. As a consequence, ω = x, which concludes the proof since s ν σ,λ (x) is odd and satisfies (reduction). We now prove (kill low SNR). Consider x > 0. Since λ → s ν σ,λ (x) is a monotonic function and s ν σ,λ (x) ∈ [0, x] for all λ, it converges for λ → 0 + to a value ω ∈ [0, x]. Assume 0 < ω < x and let 0 < ε < max(ω, x -ω). Again, we have for λ small enough

λ -ν ν,1 |t 1 | ν + λ -ν ν,2 |t 2 | ν λ -ν ν,1 |t 2 | ν + λ -ν ν,2 |t 1 | ν , (105)
⇒ (λ -ν ν,1 -λ -ν ν,2 )|t 1 | ν (λ -ν ν,1 -λ -ν ν,2 )|t 2 | ν , (106)
(x -(ω + ε)) 2 2σ 2 + λ -ν ν |ω -ε| ν < (x -t) 2 2σ 2 + λ -ν ν |t| ν . (112)
Moreover, since ω = ε, we have for λ small enough

x 2 2σ 2 < (x -(ω + ε)) 2 2σ 2 + λ -ν ν |ω -ε| ν . (113)
Combining the two last inequalities shows that

(x -0) 2 2σ 2 + λ -ν ν |0| ν < (x -t) 2 2σ 2 + λ -ν ν |t| ν , (114)
which is in contradiction with the fact that t minimizes the objective. As a consequence, ω = 0, which concludes the proof since s ν σ,λ (x) is odd and satisfies (reduction).

Figure 2 .

 2 Figure 2. Illustration of EPLL framework for image denoising with a GMM prior. A large collection of (overlapping) patches are first extracted. For each patch, an optimal Gaussian component is picked based on a measure of discrepancy with the given patch. This Gaussian component is next used as a prior model to denoise the given patch by linear shrinkage in the corresponding eigenspace and depending on β. All estimated patches are finally aggregated together, weighted by β and combined with the original noisy image to produce a first estimate. The procedure is repeated 5 times with increasing values of β.

Figure 4 .

 4 Figure 4. From left to right: the convolution of a Gaussian distribution with standard deviation λ = 2 with a Gaussian distribution with standard deviation σ = 1, the corresponding discrepancy function and logdiscrepancy function.

Figure 5 .

 5 Figure 5. From left to right: the convolution of a Laplacian distribution with standard deviation λ = 2 with a Gaussian distribution with standard deviation σ = 1, the corresponding discrepancy function and logdiscrepancy function.

Remark 3 .Figure 6 .

 36 Figure 6. From left to right: the convolution of a generalized Gaussian distribution with standard deviation λ = 2 with a Gaussian distribution with standard deviation σ = 1, the corresponding discrepancy function and log-discrepancy function. From top to bottom: the GGD has a shape parameter ν = 1.5 and .3, respectively.

Remark 4 .

 4 For ν = 1, Theorem 1 and Theorem 2 coincide with Theorem 3 and Theorem 4.

Figure 7 .

 7 Figure 7. Illustrations of the log-discrepancy function for various 0.3 ν 2 and SNR λ/σ.

Figure 8

 8 Figure 8. Illustrations of our approximations of ϕ ν λ and the corresponding underlying posterior distribution N (0, ν, λ) * G(0, 1) (where ν = .8 and λ = 4). The blue curves have been obtained by evaluating the convolution using numerical integration techniques for all x. The dashed curves are obtained using the proposed relu-and softplus-based approximations that have closed-form expressions.

Figure 9 .

 9 Figure 9. Lookup tables used to store the values of the parameters γ ν λ , β1, β2 and h for various .3 ν 2 and 10 -3 λ 10 3 . A regular grid of 100 values has been used for ν and a logarithmic grid of 100 values has been used for λ. This leads to a total of 10, 000 combinations for each of the four lookup tables.

Figure 10 .

 10 Figure 10. Illustrations of the shrinkage function for various 0.3 < ν 2 and SNR λ/σ.

  (a) GMM (ν = 2) (b) GGMM (.3 ν 2) (c) LMM (ν = 1) (d) HLMM (ν = .5)

Figure 12 .

 12 Figure 12. Average log-likelihood of all non-overlapping patches (with subtracted mean) of each of the 40 images of our validation subset of the testing BSDS dataset for the GMM, GGMM, LMM and HLMM. The total average over the 40 images is shown in the last column.

Figure 13 .

 13 Figure 13. (a) Close in on the image Castle from the BSDS testing dataset, (b) a noisy version degraded by additive white Gaussian noise with standard deviation σ = 20 and (c)-(f ) results of EPLL under four patch priors: GMM, GGMM, LMM and HLMM, respectively. PSNR and SSIM are given in the bottom-left corner.

Figure 14 .

 14 Figure 14. (a) Close in on the standard image Cameraman. (b) a noisy version degraded by additive white Gaussian noise with standard deviation σ = 20 and (c)-(f ) results of EPLL under four patch priors: GMM, GGMM, LMM and HLMM, respectively. PSNR and SSIM are given in the bottom-left corner.

Figure 15 .

 15 Figure 15. (a) Close in on the standard image Barbara. (b) a noisy version degraded by additive white noise with standard deviation σ = 20 and (c)-(f ) results of EPLL under four patch priors: GMM, GGMM, LMM and HLMM, respectively. PSNR and SSIM are given in the bottom-left corner.

Figure 16 .

 16 Figure[START_REF] Dominguez-Molina | A practical procedure to estimate the shape parameter in the generalized gaussian distribution[END_REF]. Evolution of performance of EPLL with a GGMM, LMM (ν = 1) and a GMM (ν = 2) under misspecification of the noise standard deviation σ. Performances are measured in terms of PSNR on the BSDS dataset corrupted by a Gaussian noise with standard deviation σ = 20. For each of the three priors, EPLL has been run assuming σ was ranging from 15 to 30.

29. 49 /Figure 17 .

 4917 Figure 17. Results obtained by GGMM-EPLL on a 128 × 128 cropped image of the BSDS testing dataset damaged by additive white Gaussian noise with σ = 20. These are obtained respectively by (a) evaluating the classification and shrinkage problem with numerical solvers (numerical integration and Halley's root-finding method), (b) approximating the classification problem only, and (c) approximating both problems. PSNR and SSIM are given in the bottom-left corner. Running time (averaged on ten runs) of the overall GGMM-EPLL are indicated on the captions: our accelerations lead to a speed-up of ×15, 000 and ×1.5 respectively.

Appendix D. Proof of Theorem 2 .

 2 Lemma D.1. Let a > 0 and b > 0. For x in the vicinity of +∞, we have1 2abx log erfc(ax + b) + e -4abx erfc(-ax + b) 2 = -2 + o (1) .(60)Proof. We have lim x→+∞ erfc(x) = 0 and lim x→+∞ erfc(-x) = 2, it follows that erfc(-ax + b) ∼ +∞ 2 , (61) e -4abx erfc(-ax + b) ∼ +∞ 2e -4abx , (62) erfc(ax + b) + e -4abx erfc(-ax + b) 2 ∼ +∞ e -4abx , (63) log erfc(ax + b) + e -4abx erfc(-ax + b) ax + b) + e -4abx erfc(-ax + b) ∼ +∞ -2 ,

  first-order Taylor's expansion for log(1 + x), we conclude the proof as

  v(x) = -∂ ∂x log p(x) and w(x) = -∂ ∂x log q(x) . (79) Assume v and w are positive continuous function and regularly oscillating, i.e.: )w(x) = +∞ , (81) then, for x → ∞, we have log +∞ -∞p(x -t)q(t) dt ∼ log q(x) . (82)Proof of Theorem 4. Using the definition of the discrepancy function,f ν 1,λ (x) = -log R G(t; 0, λ, ν) • N (x -t; 0, 1) dt = -log +∞ -∞

√ 2π exp - x 2 2 ⇔

 2 v(x) = -∂ ∂x log p(x) = 2x and w(x) = -x > 0 large enough and y > 0 small enoughy = p(x) ⇔ y = 1 x = -log(2π) -2 log y . )w(x) = ∞. Moreover, for ν < 2, we have lim x→∞ w(x) v(x) = lim x→∞ ν 2λ ν x ν-2 = 0 . (89)It follows that Lemma F.1 applies, and then for largex f ν 1,λ (x) ∼ -log q(x) that λ ν = λ Γ(1/ν) Γ(3/ν) , we conclude the proof since ϕ ν λ (x) = log f ν 1,λ (x) -γ ν λ ∼ ν log x -ν log λ ν .(91)Appendix G. Proof of Proposition 3.Proof. Starting from the definition of s ν σ,λ and using the change of variable t → σt, eq. (reduction) follows ass ν σ,λ (x) = argmin t∈R (x -t) 2 2σ 2 + λ -ν ν |t| ν = σargmin t∈R (x -σt) 2 2σ 2 + λ -ν ν |σt| ν , (92) = σargmin t∈R (x/σ -t) 2 2 + (λ ν /σ) -ν |t| ν = σs ν 1,λ/σ (x/σ) .(93)For eq. (odd), we use the change of variable t → -ts ν σ,λ (-x) = argmin t∈R (-x -t) 2 2σ 2 + λ -ν ν |t| ν = -argmin t∈R (-x + t) 2 2σ 2 + λ -ν ν |t| ν , (94) = -argmin t∈R (x -t) 2 2σ 2 + λ -ν ν |t| ν = -s ν σ,λ (x) . (95)

⇒ |t 1

 1 | ν |t 2 | ν (since λ 1 > λ 2 and ν > 0) . (107)We now prove (keep high SNR). Consider x > 0. Since λ → s ν σ,λ (x) is a monotonic function and s ν σ,λ (x) ∈ [0, x] for all λ, it converges for λ → ∞ to a value ω ∈ [0, x]. Assume 0 < ω < x and let 0 < ε < max(ω, x -ω). By definition of the limit, for λ big enough0 < ω -ε < t s ν σ,λ (x) < ω + ε . (108)It follows that x -t > x -(w + ε) > 0, and then(x -(ω + ε)) 2 2σ 2 + λ -ν ν |ω -ε| ν < (x -t) 2 2σ 2 + λ -ν ν |t| ν . (109)Moreover, since ω + ε = x, we have for λ big enoughλ -ν ν |x| ν < (x -(ω + ε)) 2 2σ 2 + λ -ν ν |ω -ε| ν . (110)Combining the two last inequalities shows that(x -x) 2 2σ 2 + λ -ν ν |x| ν < (x -t) 2 2σ 2 + λ -ν ν |t| ν ,(111)

Table 1

 1 Shrinkage function under generalized Gaussian priors

Table 2

 2 Image denoising performance comparison of EPLL algorithm with GMM and GGMM priors. PSNR and SSIM values are obtained on the BSDS test set (average over 60 images), on six standard images corrupted with 5 different levels of noise (average over 10 noise realizations), and finally an average over these 66 images. BM3D algorithm results are also included for reference purposes.

	σ	Algo.	BSDS barbara	camera man	hill	house	lena	mandrill Avg.
						PSNR
		BM3D	37.33	38.30	38.28	36.04 39.82 38.70	35.26	37.36
	5	GMM	37.25	37.60	38.07	35.93 38.81 38.49	35.22	37.26
		GGMM 37.33 37.73	38.12 35.95 38.94 38.52	35.23	37.33
		BM3D	33.06	34.95	34.10	31.88 36.69 35.90	30.58	33.15
	10	GMM	33.02	33.65	33.91	31.79 35.56 35.46	30.55	33.06
		GGMM 33.10 33.87	34.01 31.81 35.72 35.59	30.58	33.15
		BM3D	29.38	31.73	30.42	28.56 33.81 33.02	26.60	29.50
	20	GMM	29.36	29.76	30.16	28.46 32.77 32.40	26.60	29.42
		GGMM 29.43 30.02	30.24 28.48 33.03 32.59	26.64	29.50
		BM3D	26.28	27.97	27.16	25.89 30.69 29.81	23.07	26.38
	40	GMM	26.21	26.02	26.93	25.68 29.60 29.18	23.25	26.26
		GGMM 26.26 26.17	27.03 25.70 29.89 29.42	23.21	26.32
		BM3D	24.81	26.31	25.24	24.52 28.74 28.19	21.71	24.90
	60	GMM	24.57	23.95	25.10	24.21 27.53 27.28	21.57	24.61
		GGMM 24.64 24.03	25.17 24.25 27.80 27.52	21.50	24.67
						SSIM
		BM3D	.9619	.9643	.9613	.9508 .9571 .9436	.9588	.9614
	5	GMM	.9626	.9616	.9604 .9511 .9475 .9434	.9597	.9618
		GGMM .9628 .9617	.9602	.9507 .9469 .9425	.9593	.9620
		BM3D	.9115	.9410	.9286	.8821 .9215 .9155	.8983	.9117
	10	GMM	.9155	.9298	.9307	.8858 .8999 .9107	.9022	.9150
		GGMM .9154	.9313	.9309	.8839 .8992 .9112	.9007	.9149
		BM3D	.8236	.9036	.8685	.7789 .8741 .8763	.7943	.8260
	20	GMM	.8315	.8687	.8704 .7812 .8596 .8639	.8030	.8324
		GGMM .8297	.8715	.8699	.7766 .8629 .8669	.7991	.8308
		BM3D	.7074	.8196	.7954	.6599 .8276 .8143	.6184	.7118
	40	GMM	.7054	.7509	.7780	.6496 .8025 .7918	.6341	.7081
		GGMM .7018	.7526	.7842	.6430 .8112 .7995	.6192	.7048
		BM3D	.6375	.7581	.7496	.5859 .7956 .7784	.4993	.6427
	60	GMM	.6212	.6534	.7174	.5661 .7507 .7350	.5001	.6241
		GGMM .6174	.6544	.7266	.5592 .7622 .7438	.4782	.6207

Table 3

 3 Image denoising performance comparison of EPLL algorithm with different priors. PSNR values are obtained on the BSDS test set (average over 60 images), on six standard images corrupted with 3 different levels of noise (average over 10 noise realizations), and finally an average over these 66 images.

	σ	Algo.	BSDS barbara	camera man	hill	house	lena	mandrill Avg.
						PSNR
		GMM	37.25	37.60	38.07	35.93 38.81 38.49	35.22	37.26
	5	LMM HLMM 36.85 37.31	37.83 37.42	38.11 37.66	35.89 38.93 38.49 35.39 38.37 38.08	35.18 34.77	37.32 36.86
		GGMM 37.33	37.73	38.12 35.95 38.94 38.52	35.23	37.33
		GMM	29.36	29.76	30.16	28.46 32.77 32.40	26.60	29.42
	20	LMM HLMM 28.48 29.30	30.18 29.28	30.04 29.04	28.36 33.22 32.72 27.72 32.50 32.10	26.43 25.44	29.37 28.56
		GGMM 29.43	30.02	30.24 28.48 33.03 32.59	26.64	29.50
		GMM	24.57	23.95	25.10	24.21 27.53 27.28	21.57	24.61
	60	LMM HLMM 23.95 24.55	23.94 23.16	24.96 23.72	24.23 27.91 27.58 23.84 27.10 26.94	21.35 20.67	24.59 23.97
		GGMM 24.64 24.03	25.17 24.25 27.80 27.52	21.50	24.67

  Approx. discrepancy Approx. disc. & shrink.

	Total	10h 29m 15s 100%	2.54s (×15, 000) 100%	×1.5 ----→	1.63s (×23, 000) 100%
	Discrepancy	10h 29m 14s	1.44s (×26, 000)		1.44s (×26, 000)
		>99.99%	56.69%		88.34%
	Shrinkage	1.08s	1.08s (×1)		0.17s (×6.3)
		<0.001%	42.52%		10.43%
	Patch extraction and reprojection	0.02s <0.001%	0.02s (×1) 0.79%		0.02s (×1) 1.23%

  and(x 2 -t 2 ) 2 2σ 2 + λ -ν ν |t 2 | ν (x 2 -t 1 ) 2 2σ 2 + λ -ν ν |t 1 | ν . (99)Summing both inequalities lead to(x 1 -t 1 ) 2 + (x 2 -t 2 ) 2 (x 1 -t 2 ) 2 + (x 2 -t 1 ) 2 , (100) ⇒ -2x 1 t 1 -2x 2 t 2 -2x 1 t 2 -2x 2 t 1 ,(101)

In[START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], this is referred to as "conditional mixing weight".

, though we did not succeed in proving it, our numerical simulations also revealed a log-linear asymptotic behavior for x → ∞ in perfect agreement with the expression of α 2 and β 2 given in

Theorem 4.

Acknowledgments. The authors would like to thank Charles Dossal and Nicolas Papadakis for fruitful discussions.

Part of the experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/).

Appendix E. Proof of Theorem 3.

Proof. We first decompose the term exp -(x-t) 2 2 involved in the definition of f ν 1,λ and rewrite e xt using its power series

= -log 1 √ 2π