Yosra Marnissi

Yuling Zheng

Emilie Chouzenoux

Jean-Christophe Pesquet
email: jean-christophe@pesquet.eu

A Variational Bayesian Approach for Restoring Data Corrupted with Non-Gaussian Noise

Keywords: inverse problems, restoration, Variational Bayesian methods, parameter estimation, non-Gaussian noise, majorization, minimization

In this paper, a methodology is investigated for signal recovery in the presence of non-Gaussian noise. In contrast with regularized minimization approaches often adopted in the literature, in our algorithm the regularization parameter is reliably estimated from the observations. As the posterior density of the unknown parameters is analytically intractable, the estimation problem is derived in a variational Bayesian framework where the goal is to provide a good approximation to the posterior distribution in order to compute posterior mean estimates. Moreover, a majorization technique is employed to circumvent the difficulties raised by the intricate forms of the non-Gaussian likelihood and of the prior density. We demonstrate the potential of the proposed approach through comparisons with state-of-the-art techniques that are specifically tailored to signal recovery in the presence of mixed Poisson-Gaussian noise. Results show that the proposed approach is efficient and achieves performance comparable with other methods where the regularization parameter is manually tuned from an available ground truth.

Introduction

One of the most challenging tasks in signal processing is restoration, where one aims at providing an accurate estimate of the original signal from degraded observations. These degradations may arise due to various phenomena which are often unavoidable in practical situations. Undesirable blurring may be introduced by the atmosphere or may also stem from the intrinsic limitations of the acquisition system characterized by its point spread function. Furthermore, data can be perturbed by noise which can be viewed as parasite signal added to the information of interest, hence altering the extraction of this information. Noise may originate from various sources. On the one hand, sensors generally suffer from internal fluctuations referred to electrical or thermal noise. This type of noise is additive, independent of the signal of interest, and it can be modeled by a Gaussian distribution. On the other hand, it has been experimentally proven that in many situations, the signal of interest may suffer from noise with more complex 1 characteristics. In fact, many devices lead to measurements distorted by heteroscedastic noise i.e. the characteristics of the noise depend on the characteristics of the unknown signal [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF][START_REF] Tian | Analysis of temporal noise in CMOS photodiode active pixel sensor[END_REF][START_REF] James R Janesick | Photon transfer[END_REF][START_REF] Azzari | Gaussian-Cauchy mixture modeling for robust signal-dependent noise estimation[END_REF][START_REF] Liu | Practical signal-dependent noise parameter estimation from a single noisy image[END_REF][START_REF] Chakrabarti | Image restoration with signal-dependent camera noise[END_REF][START_REF] Boubchir | Undecimated wavelet-based Bayesian denoising in mixed Poisson-Gaussian noise with application on medical and biological images[END_REF]. For example, to better reflect the physical properties of optical communication, the noise remains additive Gaussian but its variance is assumed to be dependent on the unknown signal [START_REF] Moser | Capacity results of an optical intensity channel with input-dependent Gaussian noise[END_REF]. Signals can also be corrupted with multiplicative noise [START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF][START_REF] Buades | A non-local algorithm for image denoising[END_REF][START_REF] Dong | A convex variational model for restoring blurred images with multiplicative noise[END_REF][START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF] such as the speckle noise which commonly affects synthetic aperture radar (SAR), medical ultrasound and optical coherence tomography images [START_REF] Parrilli | A nonlocal SAR image denoising algorithm based on LLMMSE wavelet shrinkage[END_REF], as well as with impulsive noise [START_REF] Cai | Fast two-phase image deblurring under impulse noise[END_REF]. A mixture of Gaussian and impulsive noise has also been studied in [START_REF] Xiao | Restoration of images corrupted by mixed gaussian-impulse noise via ℓ 1ℓ 0 minimization[END_REF][START_REF] Yan | Restoration of images corrupted by impulse noise and mixed gaussian impulse noise using blind inpainting[END_REF]. Furthermore, in applications such as astronomy, medicine, and fluorescence microscopy where signals are acquired via photon counting devices, like CMOS and CCD cameras, the number of collected photons is related to some non-additive counting errors resulting in a shot noise [1-3, 17, 18]. The latter is non-additive, signal-dependent and it can be modeled by a Poisson distribution [START_REF] Salmon | Poisson noise reduction with non-local PCA[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Jeong | Frame-based Poisson image restoration using a proximal linearized alternating direction method[END_REF][START_REF] Harizanov | Epigraphical projection for solving least squares Anscombe transformed constrained optimization problems[END_REF][START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF][START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF][START_REF] Bonettini | New convergence results for the scaled gradient projection method[END_REF][START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF][START_REF] Bonettini | An alternating extragradient method for total variationbased image restoration from Poisson data[END_REF][START_REF] Altmann | Lidar waveformbased analysis of depth images constructed using sparse single-photon data[END_REF]. In this case, when the noise is assumed to be Poisson distributed, the implicit assumption is that Poisson noise dominates over all other noise kinds. Otherwise, the involved noise is a combination of Poisson and Gaussian (PG) components [START_REF] Roberts | Bayesian denoising/deblurring of Poisson-Gaussian corrupted data using complex wavelets[END_REF][START_REF] Mäkitalo | Noise parameter mismatch in variance stabilization, with an application to Poisson-Gaussian noise estimation[END_REF][START_REF] Mäkitalo | Poisson-Gaussian denoising using the exact unbiased inverse of the generalized Anscombe transformation[END_REF][START_REF] Mäkitalo | Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise[END_REF][START_REF] Luisier | Image denoising in mixed Poisson-Gaussian noise[END_REF][START_REF] Li | A reweighted ℓ 2 method for image restoration with Poisson and mixed Poisson-Gaussian noise[END_REF][START_REF] Lantéri | Restoration of astrophysical images -the case of Poisson data with additive Gaussian noise[END_REF][START_REF] Jezierska | A primal-dual proximal splitting approach for restoring data corrupted with Poisson-Gaussian noise[END_REF][START_REF] Jezierska | An EM approach for time-variant Poisson-Gaussian model parameter estimation[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF][START_REF] Kittisuwan | Medical image denoising using simple form of MMSE estimation in Poisson-Gaussian noise model[END_REF][START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF]. Most of the existing denoising methods only consider the noise as independent Gaussian, mainly because of the difficulties raised in handling other noise sources than the Gaussian one. In this paper, we focus on signal recovery beyond the standard additive independent Gaussian noise assumption.

State-of-the-art

Most of the existing strategies for such problems define the estimate as a minimizer of an appropriate cost function. The latter is composed of two terms: the so-called data fidelity term whose role is to make the solution consistent with the observation and the regularization term that incorporates prior information about the target signal so as to ensure the stability of the solution [START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structures and problems[END_REF]. Several algorithms have been proposed to tackle the problem of restoration for signals corrupted with non-Gaussian noise by using minimization approaches. For example, in [START_REF] Repetti | A penalized weighted least squares approach for restoring data corrupted with signal-dependent noise[END_REF], a method is proposed to restore signals degraded by a linear operator and corrupted with an additive Gaussian noise having a signal-dependent variance. An early work in [START_REF] Snyder | Image recovery from data acquired with a charge-coupled-device camera[END_REF] and more recent developments in [29, 34-36, 39, 45, 46] have proposed to restore signals corrupted with mixed PG noise using different approximations of the PG data fidelity term. In all these approaches, the regularization parameter allows a tradeoff to be performed between fidelity to the observations and the prior information. Too small values of this parameter may lead to noisy estimates while too large values yield oversmoothed solutions. Consequently, the problem of setting a proper value of the regularization parameter should be addressed carefully and may depend on both the properties of the observations and the statistics of the target signal. When ground truth is available, one can choose the value of the regularization parameter that gives the minimal residual error evaluated through some suitable metric. However, in real applications where no ground truth is available, the problem of selecting the regularization parameter remains an open issue especially in situations where the images are acquired under poor conditions i.e. when the noise level is very high. Among existing approaches dealing with regularization parameter estimation, the works in [START_REF] Ramani | Monte-Carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF][START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF][START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF][START_REF] Deledalle | Stein unbiased gradient estimator of the risk (SUGAR) for multiple parameter selection[END_REF][START_REF] Giryes | The projected GSURE for automatic parameter tuning in iterative shrinkage methods[END_REF][START_REF] Almeida | Parameter estimation for blind and non-blind deblurring using residual whiteness measures[END_REF][START_REF] Hansen | Exploiting residual information in the parameter choice for discrete ill-posed problems[END_REF] have to be mentioned. However, most of the mentioned methods were developed under the assumption of a Gaussian noise and their extension to the context of non-Gaussian noise is not easy. One can however mention the works in [START_REF] Mäkitalo | Poisson-Gaussian denoising using the exact unbiased inverse of the generalized Anscombe transformation[END_REF][START_REF] Luisier | Image denoising in mixed Poisson-Gaussian noise[END_REF] proposing efficient estimators in the context of denoising i.e. problems that do not involve linear degradation. Other approaches can be found in [START_REF] Bertero | A discrepancy principle for Poisson data[END_REF][START_REF] Zanni | Numerical methods for parameter estimation in Poisson data inversion[END_REF] proposing efficient estimates in the specific case of a Poisson likelihood.

To address the shortcomings of these methods, one can adopt the Bayesian framework. In particular, Bayesian estimation methods based on Markov Chain Monte Carlo (MCMC) sampling algorithms have been recently extended to inverse problems involving non-Gaussian noise [START_REF] Ying | A blind receiver with multiple antennas in impulsive noise with Gaussian mixtures through MCMC approaches[END_REF][START_REF] Altmann | Spectral unmixing of multispectral Lidar signals[END_REF][START_REF] Murphy | Joint Bayesian removal of impulse and background noise[END_REF][START_REF] Chaâri | Fast joint detection-estimation of evoked brain activity in event-related FMRI using a variational approach[END_REF]. However, despite good estimation performance that has been obtained, such methods remain computationally expensive for large scale problems. Another alternative approach is to rely on variational Bayesian approximation (VBA) [START_REF] Šmídl | The Variational Bayes Method in Signal Processing[END_REF][START_REF] Zheng | Efficient variational Bayesian approximation method based on subspace optimization[END_REF][START_REF] Babacan | Variational Bayesian super resolution[END_REF][START_REF] Chen | Variational Bayesian methods for multimedia problems[END_REF][START_REF] Drémeau | Boltzmann machine and mean-field approximation for structured sparse decompositions[END_REF]. Instead of simulating from the true posterior distribution, VBA approaches aim at approximating the intractable true posterior distribution with a tractable one from which the posterior mean can be easily computed. These methods can lead generally to a relatively low computational complexity when compared with sampling based algorithms.

Contributions

In this paper, we propose a VBA estimation approach for signals degraded by an arbitrary linear operator and corrupted with non-Gaussian noise. One of the main advantages of the proposed method is that it allows us to jointly estimate the original signal and the required regularization parameter from the observed data by providing good approximations of the Minimum Mean Square Estimator (MMSE) for the problem of interest. While using VBA, the main difficulty arising in the non-Gaussian case is that the involved likelihood and the prior density may have a complicated form and are not necessarily conjugate. To address this problem, a majorization technique is adopted providing a tractable VBA solution for non-conjugate distributions. Our approach allows us to employ a wide class of a priori distributions accounting for the possible sparsity of the target signal after some appropriate linear transformation. Moreover, it can be easily applied to several non Gaussian likelihoods that have been widely used [START_REF] Roberts | Bayesian denoising/deblurring of Poisson-Gaussian corrupted data using complex wavelets[END_REF][START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF][START_REF] Murtagh | Image restoration with noise suppression using a multiresolution support[END_REF][START_REF] Starck | Automatic noise estimation from the multiresolution support[END_REF]. In particular, experiments in the case of images corrupted by PG noise showcase the good performance of our approach compared with methods using the discrepancy principle for estimating the regularization parameter [START_REF] Bardsley | Regularization parameter selection methods for ill-posed Poisson maximum likelihood estimation[END_REF]. Moreover, we propose variants of our method leading to a significant reduction of the computational cost while maintaining a satisfactory restoration quality.

Outline

This paper is organized as follows. In Section 2, we formulate the considered signal recovery problem in the Bayesian framework and we present a short overview on the variational Bayesian principle. In Section 3, we present our proposed estimation method based on VBA. In Section 5, we provide simulation results together with comparisons with state-of-the-art methods in terms of image restoration performance and computation time. Finally, some conclusions are drawn in Section 6.

Problem statement 2.1 Bayesian formulation

In this paper, we consider a wide range of applications where the degradation model can be formulated as an inverse possibly ill-posed problem as follows:

(∀i ∈ {1, . . . , M }) , y i = D([Hx] i ), (1)
where x ∈ R N is the signal of interest, H ∈ R M ×N is the linear operator typically modeling a blur, a projection, or a combination of both degradations, [Hx] i denotes the i-th component of Hx , y = (y i ) 1 i M ∈ R M is the measured data and D is the noise model that may depend on the target data. The objective is to find an estimator x of x from H and y. The neg-loglikelihood φ of the observations reads

(∀x ∈ R N ) φ(x; y) = -ln p(y|x) = M i=1 φ i ([Hx] i ; y i ). (2)
Depending on the noise statistical model D, φ i may take various forms [START_REF] James R Janesick | Photon transfer[END_REF][START_REF] Cai | Fast two-phase image deblurring under impulse noise[END_REF][START_REF] Yan | Restoration of images corrupted by impulse noise and mixed gaussian impulse noise using blind inpainting[END_REF][START_REF] Mäkitalo | Poisson-Gaussian denoising using the exact unbiased inverse of the generalized Anscombe transformation[END_REF]. In particular, it reduces to a least squares function for additive Gaussian noise. In the Bayesian framework, we apply regularization by assigning a prior distribution to the data x to be recovered. In this paper, we adopt the following flexible expression of the prior density of x:

p(x | γ) = τ γ N 2κ exp -γ J j=1 D j x 2κ , (3)
where κ is a constant in (0, 1],

• denotes the ℓ 2 -norm and (D j )

1 j J ∈ (R S×N ) J where D = [D ⊤ 1 , . . . , D ⊤ J ]
⊤ is an operator that may be a matrix computing the discrete difference between neighboring pixels or a frame-analysis operator [START_REF] Pustelnik | Wavelet-based image deconvolution and reconstruction[END_REF]. The constant γ ∈ (0, +∞) can be viewed as a regularization parameter that plays a prominent role in the restoration process and τ ∈ (0, +∞) is a normalization constant independent of γ. Note that the form of the partition function for such a prior distribution follows from the fact that the associated potential is 2κ-homogeneous [START_REF] Pereyra | Maximum-a-posteriori estimation with unknown regularisation parameters[END_REF].

In this paper, the noise shape parameter κ is chosen to be fixed through empirical methods and we aim at estimating parameter γ together with x. To this end, we choose a Gamma prior for γ, i.e. p(γ) ∝ γ α-1 exp(-βγ) where α and β are positive constants (set in practice to small values to ensure a weakly informative prior).

Using the Bayes' rule, we can obtain the posterior distribution of the set of unknown variables Θ = (x, γ) given the vector of observations y:

p(Θ | y) ∝ p(y | x)p(x | γ)p(γ). (4)
However, this distribution has an intricate form. In particular, its normalization constant does not have a closed form expression. To cope with this problem, we resort to the variational Bayesian framework. The rationale of this work is to find a simple approximation to the true posterior distortion, leading to a tractable computation of the posterior mean estimate.

Variational Bayes principle

The variational Bayes approach has been first introduced in physics [START_REF] Parisi | Statistical Field Theory[END_REF]. The idea behind it is to approximate the posterior distribution p(Θ | y) with another distribution denoted by q(Θ) which is as close as possible to p(Θ | y), by minimizing the Kullback-Leibler divergence between them [START_REF] Šmídl | The Variational Bayes Method in Signal Processing[END_REF][START_REF] Pustelnik | Wavelet-based image deconvolution and reconstruction[END_REF][START_REF] Pereyra | A survey of stochastic simulation and optimization methods in signal processing[END_REF]:

q opt = argmin q KL q(Θ) p(Θ | y) , (5)
where

KL q(Θ) p(Θ | y) = q(Θ) ln q(Θ) p(Θ | y) dΘ. ( 6
)
This minimization becomes tractable if a suitable factorization structure of q(Θ) is assumed.

In particular, we assume that q(Θ) = R r=1 q(Θ r ). Hence, the optimal density approximation q opt (Θ r ) for each variable Θ r , is obtained by minimizing the KL divergence while holding the remaining densities for the rest of variables fixed. In this case, there exists an optimal solution to the optimization problem [START_REF] Liu | Practical signal-dependent noise parameter estimation from a single noisy image[END_REF] for each density q(Θ r ) 1 r R , given by the exponential of the expectation of the joint density with respect to the distribution of all the unknown parameters except the one of interest i.e.

(∀r ∈ {1, . . . , R}) q opt (Θ r ) ∝ exp ln p(y, Θ) i =r q opt (Θ i ) [START_REF] Boubchir | Undecimated wavelet-based Bayesian denoising in mixed Poisson-Gaussian noise with application on medical and biological images[END_REF] where

• i =r q(Θ i ) = • i =r q(Θ i )dΘ i .
Due to the implicit relations existing between q opt (Θ r ) 1 r R , an analytical expression of q opt (Θ) generally does not exist. Usually, these distributions are determined in an iterative way, by updating one of the separable components (q(Θ r )) 1 r R while fixing the others [START_REF] Šmídl | The Variational Bayes Method in Signal Processing[END_REF]. Applications of classical VBA approaches can be found in [62-64, 72, 73] while improved VBA algorithms have been proposed in [START_REF] Zheng | Efficient variational Bayesian approximation method based on subspace optimization[END_REF][START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF]. Once the approximate distributions are computed, the unknown parameters are then estimated by the means of the obtained distributions.

Proposed approach

Name φi(v; yi) φ ′ i (v; yi) βi(yi) Domain of validity Noise model Gaussian 1 2σ 2 (v -yi) 2 1 σ 2 (v -yi) 1 σ 2 yi ∈ R Gaussian Cauchy ln 1 + (v -yi) 2 σ 2 2(v -yi) σ 2 + (v -yi) 2 2 σ 2 yi ∈ R Cauchy Anscombe Transform 2 yi + 3 8 -v + 3 8 2 2 - 2 yi + 3 8 v + 3 8 3 8 -3/2 yi + 3 8 yi - 3 8 Poisson Generalized Anscombe Transform 2 yi + σ 2 + 3 8 -v + σ 2 + 3 8 2 2 - 2 yi + 3 8 + σ 2 v + 3 8 + σ 2 3 8 + σ 2 -3/2 yi + 3 8 + σ 2 yi - 3 8 -σ 2 Poisson-Gaussian Shifted Poisson (v + σ 2 ) -(yi + σ 2 ) ln(v + σ 2 ) 1 - yi + σ 2 v + σ 2 yi + σ 2 σ 4 yi -σ 2 Poisson-Gaussian Weighted least squares (yi -v) 2 2(σ 2 + v) + 1 2 ln(σ 2 + v) 1 2 - (yi + σ 2 ) 2 2(v + σ 2 ) 2 + 1 2(σ 2 + v) max (yi + σ 2 ) 2 σ 6 - 1 2σ 4 , 1 54(yi + σ 2 ) 4 yi ∈ R\{-σ 2 } Poisson-Gaussian
Table 1: Examples of differentiable likelihoods satisfying Assumption 3.1. φ ′ i denotes the first derivative of function φ i and β i (y i ) is the Lipschitz constant of φ ′ i (for functions in lines 3-6, we assume that φ i is replaced on R -by its quadratic extension [START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF].) The expression for the Lipschitz constant of the gradient of the weighted least squares likelihood was established in [START_REF] Repetti | Algorithmes d'optimisation en grande dimension: applications à la résolution de problèmes inverses[END_REF]Chap. IV].

In this work, we assume the following separable form for q: q(Θ) = q(x)q(γ).

Unfortunately, by using directly [START_REF] Boubchir | Undecimated wavelet-based Bayesian denoising in mixed Poisson-Gaussian noise with application on medical and biological images[END_REF], we cannot obtain an explicit expression of q(x) due to the intricate form of both the prior distribution and the likelihood when the statistics of the noise are no longer Gaussian. In this paper, we propose to use deterministic methods to construct quadratic upper bounds for the negative logarithms of both the likelihood and the prior density [START_REF] Seeger | Fast variational Bayesian inference for non-conjugate matrix factorization models[END_REF]. This allows us to derive an upper bound of the desired cost function in (5) as will be described in the following.

Construction of the majorizing approximation

Likelihood

One popular approach in signal recovery is the half-quadratic formulation [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF]. Under some mild assumptions and by introducing some auxiliary variables, a complicated criterion can be written as the infimum of a surrogate half-quadratic function i.e. the latter is quadratic with respect to the original variables and the auxiliary variables appear decoupled. This half-quadratic criterion can be then efficiently minimized using classical optimization algorithms. This formulation has been widely used in energy-minimization approaches [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF][START_REF] Nikolova | Analysis of half-quadratic minimization methods for signal and image recovery[END_REF][START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF] where the initial optimization problem is replaced by the minimization of the constructed surrogate function. Furthermore, this technique has been recently extended to sampling algorithms [START_REF] Marnissi | An auxiliary variable method for langevin based MCMC algorithms[END_REF]. The initial intractable posterior distribution to sample from is replaced by the conditional distribution of the target signal given the auxiliary variables. The obtained distribution has been shown to be much simpler to explore by using standard sampling algorithms. In this paper, we propose to use half-quadratic approaches to construct an upper bound for the objective function in [START_REF] Liu | Practical signal-dependent noise parameter estimation from a single noisy image[END_REF].

We assume that the likelihood satisfies the following property:

Assumption 3.1. For every i ∈ {1, . . . , M }, φ i is differentiable on R and there exists µ i (y i ) > 0 such that the function defined by

v → v 2 2 -φ i (v;y i ) µ i (y i ) is convex on R.
In particular, this assumption is satisfied when, for every i ∈ {1, . . . , M }, φ i is β i (y i )-Lipschitz differentiable on R for every µ i (y i ) β i (y i ). Table 1 shows some examples of useful functions satisfying the desired property (up to an additive constant). Note that, since the functions in lines 3-6 of Table 1 are β i (y i )-Lipschitz differentiable only on R + , we propose to use on R -a quadratic extension of them defined as follows:

(∀v ∈ R -) φ i (v; y i ) = φ i (0; y i ) + φ ′ i (0; y i )v + 1 2 φ ′′ i (0; y i )v 2 , (9)
so that the extended version of φ i (.;

y i ) is now differentiable on R with β i (y i )-Lipschitzian gradient.
For every i ∈ {1, . . . , M } and v ∈ R, let us define the following function:

(∀v ∈ R) ς i (v; y i ) = sup t∈R - 1 2 (v -t) 2 + φ i (t; y i ) µ i (y i ) . ( 10
)
Then, the following property holds:

Proposition 3.1. For every i ∈ {1, . . . , M },

(∀v ∈ R) φ i (v; y i ) = inf w i ∈R T i (v, w i ; y i ). ( 11
)
where, for every v ∈ R,

T i (v, w i ; y i ) = µ i (y i ) 1 2 (v -w i ) 2 + ς i (w i ; y i ) . (12)
Moreover, the unique minimizer of

w i → T i (v, w i ; y i ) reads w i (v) = v - 1 µ i (y i ) φ ′ i (v; y i ). ( 13
)
Proof. See Appendix 6.

It follows from this result that

∀x ∈ R N φ(x; y) = inf w∈R M T (x, w; y), (14)
where T (x, w; y)

= M i=1 T i ([Hx] i , w i ; y i ).
Note that [START_REF] Dong | A convex variational model for restoring blurred images with multiplicative noise[END_REF] shows that, for every i ∈ {1, . . . , M }, φ i (•; y i ) is a so-called Moreau envelope of the function ς i (•; y i ). A more direct proof of Proposition 3.1 can thus be derived from the properties of the proximity operator [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] when the functions (φ i ) 1≤i≤M are convex. The proof we provide in the appendix however does not make such a restrictive assumption.

Prior

Similarly, we construct a surrogate function for the prior distribution. More precisely, we follow the same idea as in [START_REF] Chen | Variational Bayesian methods for multimedia problems[END_REF] and we use the following convexity inequality to derive a majorant for the ℓ κ -norm :

(∀ν ≥ 0)(∀υ > 0) υ κ ≤ (1 -κ)ν κ + κν κ-1 υ.
Hence, we obtain the following majorant function for the negative logarithm of the prior distribution:

γ J j=1 D j x 2κ γ J j=1 κ D j x 2 + (1 -κ)λ j λ 1-κ j . ( 15
)
where (λ j ) 1 j J are positive variables. In the following, we will denote by Q(x, λ; γ)

= J j=1 Q j (D j x, λ j ; γ),
the function in the right-hand side of the above inequality where, for every j ∈ {1, . . . , J},

Q j (D j x, λ j ; γ) = γ κ D j x 2 + (1 -κ)λ j λ 1-κ j . (16)

Proposed majorant

Thus, we can derive the following lower bound for the posterior distribution:

p(Θ | y) L(Θ|y; w, λ), (17)
where function L is defined as

L(Θ|y; w, λ) = C(y) exp [-T (x, w; y) -Q(x, λ; γ)] p(γ) with C(y) = p(y) -1 (2π) -M/2 τ γ N 2κ
. The minorization of the distribution leads to an upper bound for the KL divergence:

KL(q(Θ) p(Θ | y)) KL(q(Θ) L(Θ|y; w, λ)). ( 18
)
Note that, although the constructed lower bound in ( 17) is tangent to the posterior distribution i.e.

p(Θ | y) = sup

w∈R M ,λ∈R J L(Θ|y; w, λ), the tangency property may not be generally satisfied in [START_REF] Petropulu | Power-law shot noise and its relationship to long-memory α-stable processes[END_REF]. Thus, the tightness of the constructed majorant of the KL divergence may have a significant impact on the accuracy of the method. By minimizing the constructed bound [START_REF] Petropulu | Power-law shot noise and its relationship to long-memory α-stable processes[END_REF] with respect to w and λ, we make this bound as tight as possible. Note that, for every i ∈ {1, . . . , M } and j ∈ {1, . . . , J}, λ j → KL(q(Θ) L(Θ|y; w, λ)) and w i → KL(q(Θ) L(Θ|y; w, λ)) can be minimized separately. Hence, Problem (5) can be solved by the following four-step alternating optimization scheme:

• Minimizing the upper bound in (18) w.r.t. q(x);

• Updating the auxiliary variables w i in order to minimize KL(q(Θ) L(Θ|y; w, λ)), for every i ∈ {1, . . . , M };

• Updating the auxiliary variable λ j in order to minimize KL(q(Θ) L(Θ|y; w, λ)), for every j ∈ {1, . . . , J};

• Mimimizing the upper bound in (18) w.r.t. q(γ).

The main benefit of this majorization strategy is to guarantee that the optimal approximate posterior distribution for x belongs to the Gaussian family and the optimal approximate posterior distribution for γ belongs to the Gamma one, i.e.

q(x) ≡ N (m, Σ), q(γ) ≡ G(a, b).
Therefore, the distribution updates can be performed by updating their parameters, namely m, Σ, a, and b.

Iterative algorithm

Subsequently, at a given iteration k of the proposed algorithm, the corresponding estimated variables will be indexed by k.

Updating q(x)

Because of the majorization step, we need to minimize the upper bound on the KL divergence.

The standard solution [START_REF] Boubchir | Undecimated wavelet-based Bayesian denoising in mixed Poisson-Gaussian noise with application on medical and biological images[END_REF] can still be used by replacing the joint distribution by a lower bound L(Θ, y; w, λ) chosen proportional to L(Θ|y; w, λ):

q k+1 (x) ∝ exp ln L(x, γ, y; w k , λ k ) q k (γ) ∝ exp ln L(x, γ, y; w k , λ k )q k (γ)dγ ∝ exp - M i=1 1 2 µ i (y i ) [Hx] i -w k i 2 - a k b k J j=1 κ D j x 2 + (1 -κ)λ k j (λ k j ) 1-κ . ( 19
)
The above distribution can be identified as a multivariate Gaussian distribution whose covariance matrix and mean parameter are given by

Σ -1 k+1 = H ⊤ Diag(µ(y))H + 2 a k b k D ⊤ Λ k D, (20)
m k+1 = Σ k+1 H ⊤ u, (21)
where µ(y) = [µ 1 (y 1 ), . . . , µ M (y M )] ⊤ , u is a M × 1 vector whose i-th component is given by u i = µ i (y i )w k i and Λ is the diagonal matrix whose diagonal elements are κ(λ k j ) κ-1 I S 1≤j≤J .

Updating w

The auxiliary variable w is determined by minimizing the upper bound of KL divergence with respect to this variable:

w k+1 = argmin w q k+1 (x)q k (γ) ln q k+1 (x)q k (γ) L(Θ|y; w, λ k ) dxdγ = argmin w q k+1 (x)q k (γ) -ln L(Θ|y; w, λ k ) dxdγ = argmin w q k+1 (x) M i=1 T i ([Hx] i , w i ; y i )dx (22)
= argmin

w M i=1 T i ([Hm k+1 ] i , w i ; y i ), (23)
where the equality in [START_REF] Harizanov | Epigraphical projection for solving least squares Anscombe transformed constrained optimization problems[END_REF] follows from the expression in [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF]. Interestingly, it follows from Property 3.1 that

w k+1 i = argmin w i T i ([Hm k+1 ] i , w i ; y i ) = [Hm k+1 ] i - 1 µ i (y i ) φ ′ i ([Hm k+1 ] i ; y i ). ( 24
)

Updating λ

The variable λ is determined in a similar way: for every j ∈ {1, . . . , J}, λ k+1 j = argmin

λ j ∈[0,+∞)
KL(q k+1 x)q k (γ) L(Θ|y; w k+1 , λ)

= argmin

λ j ∈[0,+∞) Q i=1 q k+1 (x)q k (γ)Q i (D i x, λ i ; γ)dxdγ = argmin λ j ∈[0,+∞) q k+1 (x)q k (γ)Q j (D j x, λ j ; γ)dxdγ = argmin λ j ∈[0,+∞) q k+1 (x)q k (γ)γ κ D j x 2 + (1 -κ)λ j λ 1-κ j dxdγ = argmin λ j ∈[0,+∞) κ E q k+1 (x) D j x 2 + (1 -κ)λ j λ 1-κ j . ( 25
)
The minimum is attained at

λ k+1 j = E q k+1 (x) D j x 2 = D j m k+1 2 + trace D ⊤ j D j Σ k+1 . (26)
4.4 Updating q(γ)

Using [START_REF] Boubchir | Undecimated wavelet-based Bayesian denoising in mixed Poisson-Gaussian noise with application on medical and biological images[END_REF] where the joint distribution is replaced by its lower bound function, we obtain q k+1 (γ) ∝ exp ln L(x, γ, y; w k+1 , λ k+1 )

q k+1 (x)
∝ exp ln L(x, γ, y; w k+1 , λ k+1 )q k+1 (x)dx

∝γ N 2κ +α-1 exp(-βγ) exp -γ J j=1 κE q k+1 (x) D j x 2 + (1 -κ)λ k+1 j (λ k+1 j ) 1-κ ≡G(a k+1 , b k+1 ). ( 27
)
Using [START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF], one can recognize that the above distribution is a Gamma one with parameters

a k+1 = N 2κ + α = a, b k+1 = J j=1 (λ k+1 j ) κ + β. (28)

Resulting algorithm

The proposed method is outlined in Algorithm 2. It alternates between the update of the auxiliary variables and the distribution of the unknown parameters.

Algorithm 1 VBA approach for recovery of signals corrupted with non-Gaussian noise.

1. Set initial values: w 0 , λ 0 , b 0 . Compute a with (28). 2. For k = 0, 1, . . .

(a) Update parameters Σ k+1 and m k+1 of q k+1 (x) using ( 20) and ( 21). (b) Update w k+1 using [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. (c) Update λ k+1 using ( 26). (d) Update parameter b k+1 of q k+1 (γ) using (28).

Implementation issues

An additional difficulty arising in the implementation of Algorithm 2 is that the determination of Σ k+1 requires inverting the matrix given by [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF], which is computationally expensive in high dimension. To bypass this operation, we propose to compare two approaches. The first one follows the idea in [START_REF] Babacan | Variational Bayesian super resolution[END_REF]: we make use of the linear conjugate gradient method to approximate m k+1 iteratively and in [START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF], where an explicit form of Σ k+1 cannot be sidestepped, this matrix is approximated by a diagonal one whose diagonal entries are equal to the inverse of the diagonal elements of Σ -1 k+1 . The second technique uses Monte-Carlo sample averaging to approximate m k+1 and λ k+1 j : specifically, we generate samples (n s ) 1≤s≤N s from Gaussian distribution with mean m k+1 and covariance matrix Σ k+1 using [START_REF] Orieux | Sampling high-dimensional Gaussian distributions for general linear inverse problems[END_REF], as summarized in Algorithm 2. This estimator has two desirable properties. First, its accuracy is independent of the problem size, its relative error only depends on the number of samples and it decays as 2/N s (only N s = 2/ρ 2 samples are required to reach a desired relative error ρ) [START_REF] Papandreou | Gaussian sampling by local perturbations[END_REF]. Second, for the simulation of N s independent Gaussian samples, one can take advantage of a multiprocessor architecture by resorting to parallel implementation allowing us to reduce the computation time.

Algorithm 2 Stochastic approach for computing the parameters of q(x). 1. For s = 1, 2, . . . , N s (a) Perturbation : Generate

ν s ∼ N u, Diag µ(y) 1/2 η s ∼ N 0, 2γ k Λ 1/2 k with γ k = a k /b k . (b) Optimization:
Compute n s as the minimizer of

J (v) = ν s - Diag µ(y) Hv 2 (Diag(µ(y))) -1 + 1 2γ k η s -2γ k Λ k Dv 2 Λ -1 k , which is equivalent to minimize J (v) = v ⊤ Σ -1 k+1 v -2v ⊤ z s where z s = H ⊤ ν s + D ⊤ η s .
The minimizer is computed using the conjugate gradient algorithm.

Update

m k+1 = 1 N s N s s=1 n s (∀j ∈ {1, . . . , J}) λ k+1 j = 1 N s N s s=1 D j n s 2 .
5 Application to Poisson-Gaussian image restoration

Problem formulation

In this section, we propose to illustrate the usefulness of our algorithm via experiments on PG noise models in the context of image restoration. More specifically, the vector of observations y = (y i ) 1≤i≤M ∈ R M is related to the original image x through

y = z + b, (29)
where z and b are assumed to be mutually independent random vectors and

z | x ∼ P(Hx), b ∼ N (0, σ 2 I M ),
P denoting the independent Poisson distribution. Hence, we can deduce the following expression of the likelihood function [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF]:

p(y | x) = M i=1 +∞ n=1 e -[Hx] i ([Hx] i ) n n! e -1 2σ 2 (y i -n) 2 √ 2πσ 2 . ( 30
)
The exact expression of the PG likelihood (30) involves an infinite sum which makes its exact computation impossible in practice. In [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF], the infinite sum was replaced by a finite summation with bounds depending on the current estimate of x. However, this method implies a higher computational burden in the reconstruction process when compared with other likelihoods proposed in the literature as accurate approximations of [START_REF] Mäkitalo | Noise parameter mismatch in variance stabilization, with an application to Poisson-Gaussian noise estimation[END_REF]. In [START_REF] Marnissi | Fast variational Bayesian signal recovery in the presence of Poisson-Gaussian noise[END_REF], VBA inference techniques have been successfully applied to the restoration of data corrupted with PG noise using the generalized Anscombe transform (GAST) likelihood [START_REF] Mäkitalo | Poisson-Gaussian denoising using the exact unbiased inverse of the generalized Anscombe transformation[END_REF][START_REF] Mäkitalo | Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise[END_REF][START_REF] Murtagh | Image restoration with noise suppression using a multiresolution support[END_REF][START_REF] Starck | Automatic noise estimation from the multiresolution support[END_REF]. In this paper, we will employ the GAST likelihood , as well as the shifted Poisson likelihood (SPoiss) [START_REF] Chakrabarti | Image restoration with signal-dependent camera noise[END_REF] and the weighted least squares approximation (WL2) [START_REF] Li | A reweighted ℓ 2 method for image restoration with Poisson and mixed Poisson-Gaussian noise[END_REF][START_REF] Repetti | A penalized weighted least squares approach for restoring data corrupted with signal-dependent noise[END_REF][START_REF] Benvenuto | The study of an iterative method for the reconstruction of images corrupted by Poisson and Gaussian noise[END_REF] (see Table 1).

Numerical results

We evaluate the performance of the proposed approach for the restoration of images degraded by both blur and PG noise. We consider six test images, displayed in Figure 1, whose intensities have been rescaled so that pixel values belong to a chosen interval [0, x + ]. Images x 1 and x 6 are HST astronomical images while images x 2 , x 3 , x 4 and x 5 correspond to the set of confocal microscopy images considered in [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF]. These images are then artificially degraded by an operator H modeling spatially invariant blur with point spread function h and by PG noise with variance σ 2 . In our experiments, we choose a total variation prior, i.e. κ = 1/2 and for every pixel j ∈ {1, . . . , N },

D j x = [∇ h x] j , [∇ v x] j ⊤ ∈ R 2
where ∇ h and ∇ v are the discrete gradients computed in the horizontal and vertical directions. As a result, J = N . It can be noted that we can apply the proposed algorithm to all approximations of the PG fidelity term presented in Table 1. For the GAST and the SPoiss approximations, we use µ i (y i ) ≡ max {β i (y i ), ε} where ε > 0. For the WL2 approximation, we use µ i (y i ) = max (y i + σ 2 ) 2 /σ 6 , ε , so as to satisfy Assumption 3.1. Note that a data truncation is performed as a pre-processing step on the observed image y in order to satisfy the domain condition given in the fifth column of Table 1.

The goal of our experiments is twofold. First, for each likelihood, we compare the accuracy of the two proposed approximations of the covariance matrix described in Section 4.6 namely the diagonal approximation (denoted as approximation 1) and the Monte Carlo averaging strategy (denoted as approximation 2) with different number of samples N s , namely N s = 160 or 640. Second, the proposed method is compared with state-of-the-art algorithms that compute the MAP estimate. More specifically, as GAST and SPoiss likelihoods are convex and Lipschitz differentiable, we use the method presented in [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF] where a primal-dual splitting algorithm is adopted to minimize the associated penalized criterion. For the WL2 approximation, we consider the variable metric forward-backward algorithm presented in [START_REF] Repetti | A penalized weighted least squares approach for restoring data corrupted with signal-dependent noise[END_REF]. For all these approaches, we test two variants. First, we propose to estimate the regularization parameter using the discrepancy principle [START_REF] Bertero | A discrepancy principle for Poisson data[END_REF][START_REF] Zanni | Numerical methods for parameter estimation in Poisson data inversion[END_REF][START_REF] Bardsley | Regularization parameter selection methods for ill-posed Poisson maximum likelihood estimation[END_REF]. In the second variant, the regularization parameter γ is adjusted empirically to achieve the maximum Signal-to-Noise Ratio (SNR) value defined as SNR= 20 log 10 x x-x , which requires the availability of the true image. Tables 234567report the results obtained with the different images in terms of SNR, SSIM [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], and approximate computation time needed for convergence. For each likelihood, we emphasize in bold the approximation of the covariance matrix that achieves the best quantitative result in the shortest computational time. Simulations are performed on an Intel(R) Xeon(R) CPU E5-2630, @ 2.40 GHz, using a Matlab7 implementation. All tested methods are initialized with the degraded image. The Monte Carlo averaging approximation is computed using parallel implementation with 16 cores by means of the command PARFOR of the Matlab R Parallel Computing Toolbox TM .

One can observe that in most studied situations (see Tables 4567), the diagonal approximation of the covariance matrix appears to give satisfactory qualitative results in a small computation time. However, in few other situations (see Tables 2 and3), it fails to capture the real qualitative structures of the covariance matrix leading to a poorer performance. The latter issue is well alleviated by using the Monte Carlo approximation where good results, in terms of image quality, are achieved within N s = 160 samples which is equivalent to a relative approximation error equal to 11%. A few improvements are observed by decreasing the approximation error to 5% using N s = 640 samples.

We also notice that the GAST approximation does not seem to be suitable for very low count images (see Tables 2 and3), whereas, the other likelihoods lead to competitive results in all the experiments. The best tradeoff between restoration quality and small computational time seems to be achieved by the WL2 approximation. Finally, it can be observed that, in Tables 2 and4, our VBA method yields comparable performance in terms of SNR to the MAP estimate when the latter is computed with the optimal regularization parameter, while our approach requires less time to converge. In the other experiments, our approach leads to the best qualitative results. For instance, in Table 3, the gain in terms of SNR reaches up to 0.20 dB compared with the MAP estimator using the best regularization parameter, but our approach needs more time to converge. In Tables 5,6 and 7, we achieve both the best quantitative results and the smallest computational time. It should be noted that for most tested scenarii, discrepancy based approaches perform relatively poorly compared with the other methods, especially in the case of low count images (see Table 2).

In Figures 234567, we show some examples of visual results obtained with the different approaches, when the best approximation strategy for the covariance matrix is retained in the VBA method. It can be noticed that, unlike the other methods, the reconstructed images with the proposed VBA algorithm exhibit very few artifacts, without over-smoothed regions.

It should be emphasized that the problem of setting the regularization parameter for these algorithms must be carefully addressed as it highly impacts the quality of the restored image. The main advantage of our approach is that this parameter is tuned automatically without the need of the ground truth, while also often being the most competitive in terms of computation time. Furthermore, the performance of the proposed method could be further improved by using parallel implementation with more than 16 cores for the Monte Carlo approximation of the covariance matrix allowing either generating a higher number of samples (i.e. an improved estimation error) or a reduction of the computation time.

Conclusion

In this paper, we have proposed a variational Bayesian approach for solving signal recovery problems in the presence of non-Gaussian noise. Our approach has two main advantages. First, the regularization parameter is tuned automatically during the recovery process. Second, the designed method is applicable to a wide range of prior distributions and data fidelity terms. Simulations carried out on various images corrupted with mixed Poisson-Gaussian noise have shown that the proposed strategy constitutes a competitive solution for low computational and high-quality restoration of images compared with state-of-the art methods.

[Proof of Proposition 3.1] Let i ∈ {1, . . . , M }. Let us define g i : R → R such that

(∀v ∈ R) g i (v) = v 2 2 - φ i (v; y i ) µ i (y i ) . ( 31
)
According to Assumption 3.1, g i is convex. Its conjugate function [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Chapter 13] reads:

(∀w ∈ R) g * i (w) = sup v∈R (vw -g i (v)) (32)
= sup

v∈R vw + φ i (v; y i ) µ i (y i ) - v 2 2 (33) = sup v∈R - 1 2 (v -w) 2 + φ i (v; y i ) µ i (y i ) + w 2 2 . ( 34
)
According to Definition [START_REF] Buades | A non-local algorithm for image denoising[END_REF],

(∀w ∈ R) g * i (w) = ς i (w; y i ) + w 2 2 . ( 35
)
The conjugate of g * i is

(∀v ∈ R) g * * i (v) = sup w∈R vw -g * i (w) (36)
= sup

w∈R vw - w 2 2 -ς i (w; y i ) = sup w∈R - 1 2 (v -w) 2 -ς i (w; y i ) + v 2 2 = -inf w∈R 1 2 (v -w) 2 + ς i (w; y i ) + v 2 2 . ( 37
)
Since g i is convex,

g i = g * * i so that (∀v ∈ R) - φ i (v; y i ) µ i (y i ) = -inf w∈R 1 2 (v -w) 2 + ς i (w; y i ) (38)
which is equivalent to

(∀v ∈ R) φ i (v; y i ) = µ i (y i ) inf w∈R 1 2 (v -w) 2 + ς i (w; y i ) , (39)
so that (11) holds.

For every v ∈ R, let

w i (v) = g ′ i (v). (40)
By using standard properties of the subdifferential of convex functions [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Corollary 16.24], the latter equality is equivalent to

v ∈ ∂g * i w i (v) , (41)
which is a necessary and sufficient condition for w i (v) to be a minimizer of the convex function w → g * i (w)vw. This minimizer is unique since w i (v) is uniquely defined by [START_REF] Kittisuwan | Medical image denoising using simple form of MMSE estimation in Poisson-Gaussian noise model[END_REF]. We have therefore established that

w i (v) = argmax w∈R vw -g * i (w) . (42)
The definition of g i in [START_REF] Mäkitalo | Poisson-Gaussian denoising using the exact unbiased inverse of the generalized Anscombe transformation[END_REF] shows that (40) also reads

w i (v) = v - 1 µ i (y i ) φ ′ i (v; y i ). ( 43
)
According to [START_REF] Jezierska | A primal-dual proximal splitting approach for restoring data corrupted with Poisson-Gaussian noise[END_REF], it is straightforward that w i (v) also reaches the infimum in [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF]. Hence the result by using ( 43) and ( 12).

Figure 1 :

 1 Figure 1: Original images.

( a )

 a Degraded image with SNR= -2.55 dB (Uniform kernel 5 × 5, x + = 10 and σ 2 = 4). (b) Restored image with VBA approach using the Monte Carlo approximation with 640 samples: SNR= 10.27 dB (c) Restored image with discrepancy principle: SNR= 10.17 dB (d) Restored image with best parameter: SNR= 10.39 dB

Figure 2 :

 2 Figure 2: Restoration results for image x1 using WL2 approximation.

( a )

 a Degraded image with SNR= 2.21 dB (Gaussian kernel 25×25, std 1.6, x + = 12 and σ 2 = 9). (b) Restored image with VBA approach using the Monte Carlo approximation with 640 samples: SNR= 19.12 dB (c) Restored image with discrepancy principle: SNR= 17.41 dB (d) Restored image with best parameter: SNR= 18.73 dB

Figure 3 :

 3 Figure 3: Restoration results for image x2 using SPoiss approximation.

( a )

 a Degraded image with SNR= 3.14 dB (Uniform kernel 5 × 5, x + = 15 and σ 2 = 9). (b) Restored image with VBA approach using the Monte Carlo approximation with 640 samples: SNR= 12.36 dB (c) Restored image with discrepancy principle: SNR= 12.38 dB (d) Restored image with best parameter: SNR= 12.45 dB

Figure 4 :

 4 Figure 4: Restoration results for image x3 using SPoiss approximation.

( a )

 a Degraded image with SNR=7.64 dB (Uniform kernel 5 × 5, x + = 20 and σ 2 = 9). (b) Restored image with VBA approach using the diagonal approximation: SNR= 13.80 dB (c) Restored image with discrepancy principle: SNR= 13.48 dB (d) Restored image with best parameter: SNR= 13.60 dB

Figure 5 :

 5 Figure 5: Restoration results for image x4 using GAST approximation.

Table 2 :

 2 Restoration results for image x 1 with x + = 10 and σ 2 = 4. Uniform kernel with size 5 × 5. Initial SNR= -2.55 dB.

	GAST	SPoiss	WL2

Table 3 :

 3 Restoration results for the image x 2 , x + = 12 and σ 2 = 9. Gaussian kernel with size 25 × 25, std 1.6. Initial SNR= 2.21 dB.

	GAST	SPoiss	WL2

Table 4 :

 4 Restoration results for the image x 3 with x + = 15 and σ 2 = 9. Uniform kernel with size 5 × 5. Initial SNR= 3.14 dB.

	GAST	SPoiss	WL2

Table 5 :

 5 Restoration results for the image x 4 with x + = 20 and σ 2 = 9. Uniform kernel with size 5 × 5. Initial SNR= 7.64 dB.

	GAST	SPoiss	WL2

Table 6 :

 6 Restoration results for image x 5 with x + = 20 and σ 2 = 9. Gaussian kernel with size 7 × 7, std 1. Initial SNR= 8.55 dB.

	GAST	SPoiss	WL2

Table 7 :

 7 Restoration results for the image x 6 with x + = 100 and σ 2 = 36. Uniform kernel with size 3 × 3. Initial SNR= 10.68 dB.

(a) Degraded image with SNR= 8.55 dB (Gaussian kernel 7 × 7, std 1. x + = 20 and σ 2 = 9). (b) Restored image with VBA approach using the diagonal approximation: SNR= 20.71 dB (c) Restored image with discrepancy principle: SNR= 18.70 dB (d) Restored image with best parameter: SNR= 20.44 dBFigure 6: Restoration results for image x5 using WL2 approximation.

This work was supported by the CNRS Imag'in project under grant 2015OPTIMISME, by the CNRS MASTODONS project under grant 2016TABASCO, and by the Agence Nationale de la Recherche under grant ANR-14-CE27-0001 GRAPHSIP.