A Fraysse

T Rodet
email: trodet@satie.ens.cachan.fr

A MEASURE-THEORETIC VARIATIONAL BAYESIAN ALGORITHM FOR LARGE DIMENSIONAL PROBLEMS

Keywords: ill-posed inverse problems, variational Bayesian methodology, sparse signal reconstruction, infinite dimensional convex optimization

In this paper we provide an algorithm adapted to the variational Bayesian approximation. The main contribution is to transpose a classical iterative algorithm of optimization in the metric space of measures involved in the Bayesian methodology. Once given the convergence properties of this algorithm, we consider its application to large dimensional inverse problems, especially for unsupervised reconstruction. The interest of our method is enhanced by its application to large dimensional linear inverse problems involving sparse objects. Finally, we provide simulation results. First we show the good numerical performances of our method compared to classical ones on a small example. On a second time we deal with a large dimensional dictionary learning problem.

1. Introduction. The recent development of information technologies has increased the expansion of inverse problems for very large dimensional datasets. There is currently a significant growth in the number of measurements involved in reconstruction problems together with an expansion in the size of considered objects. Thus it is often necessary to handle the reconstruction of more than one million parameters. At the same time, signal processing techniques have helped to overcome the limitations of measurement instruments and to supply the design of systems involving indirect measures. These new instruments introduced different signal processing challenges, such as super resolution, deconvolution, source separation or tomographic reconstruction. All these problems are ill posed, the only information contained in the data and in the model of acquisition are not sufficient to obtain a good estimation of the source objects.

To solve these ill-posed problems, additional a priori information is often needed. The Bayesian approach appearing in this paper consists in a modelisation of sources of information as probability density functions [START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structure and problems[END_REF][START_REF] Mackay | Information theory, inference, and learning algorithms[END_REF][START_REF]Bayesian approach to inverse problems[END_REF]. This approach allows the development of unsupervised estimation methods, such that the parameters of probability distributions also called hyperparameters, typically the different variances, are adjusted automatically. These hyperparameters can tune the trade-off between information coming from data (likelihood) and a priori information. All the information related to the object and to the hyperparameters is summarized by the posterior distribution. This distribution allows the construction of several estimators such as the maximum a posteriori (MAP) or the posterior mean (PM), optimal for some given cost functions. These estimators can in fact be obtained by summing up the information contained in the posterior distribution. However the information given only by these estimators cannot be sufficient in the context of inverse problems, especially when recursive or online methods are required. Moreover the determination of these estimators needs an exact knowledge of the posterior distribution, which is generally intractable. Hence this posterior distribution is in practice approximated thanks, in general, to numerical methods. The most classical method is given by Markov Chains Monte Carlo(MCMC) methods, see [START_REF] Robert | Simulations par la méthode MCMC[END_REF][START_REF] Robert | Monte-Carlo statistical methods[END_REF], or particular filtering [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF] in recursive frameworks. However, these methods are no longer tractable for large dimensional problems. Indeed, in this case, the rate of convergence is too slow to allow reasonable practical use.

Another approach, considered here, is to determine an analytical approximation which leads to Laplace approximation or to Variational Bayesian Approximation (VBA), [START_REF]Variational bayesian filtering[END_REF]. These approaches are more suitable when the full information has to be propagated along recursive methods, see [START_REF]Variational bayesian filtering[END_REF][START_REF] El Fquih | Variational bayesian Kalman filtering in dynamical tomography[END_REF] or for experimental design [START_REF] Seeger | Bayesian inference and optimal design for sparse linear model[END_REF][START_REF] Seeger | Large scale bayesian inference and experimental design for sparse linear models[END_REF] where marginal distributions need to be calculated quickly. One of the first works on variational Bayesian inference was introduced by D. MacKay in [START_REF] Mackay | Ensemble learning and evidence maximization[END_REF]. This methodology has been used in a computer science context, see [START_REF] Jordan | An Introduction to variational Methods for Graphical Models[END_REF] for instance, and in signal processing for different applications such as source separation using ICA [START_REF] Miskin | Ensemble learning for independent component analysis[END_REF][START_REF] Choudrey | Variational methods for bayesian independent component analysis[END_REF], Mixture Models estimation [START_REF] Morton | Variational bayesian learning for mixture autoregressive models with uncertain-order[END_REF], hyperspectral imaging reduction [START_REF] Bali | Bayesian approach with hidden Markov modeling and mean field approximation for hyperspectral data analysis[END_REF], deconvolution [START_REF] Chantas | Variational Bayesian image restoration based on a product of t-distributions image prior[END_REF][START_REF] Babacan | Variationnal Bayesian Blind Deconvolution Using a Total Variation Prior[END_REF] or recursive methods [START_REF] Smidl | The variational bayes method in signal processing[END_REF][START_REF]Variational bayesian filtering[END_REF]. The main outline of the variational Bayesian methodology is to determine an approximation, in terms of Kullback-Leiber divergence, of the posterior distribution by a simpler, typically separable, law. Thanks to this method, the initial inverse problem appears as a convex optimization problem in a function space. Moreover, in classical Bayesian variational methods, the prior distribution is chosen to be conjugated with the likelihood one. Hence the posterior distribution belongs to a given family and the optimization becomes an optimization of its parameters. Unfortunately, this variational method induces intricate formula that require the implementation of numerical approximations. They are hence approximated by fixed point methods such as alternate minimization. For large dimensional problems this turn out to be an important drawback. Some recent works, such as [START_REF]Variational bayesian super resolution[END_REF] [39], proposed accelerated variational Bayesian algorithms based on iterative approximations for update equations.

In the present paper we choose a different approach by directly considering a numerical resolution of the infinite dimensional optimization problem involved in the variational Bayesian framework. This approach induces fewer approximations and ensures numerical convergence to the true solution of the variational Bayesian optimization problem. The goal is to construct an iterative algorithm able to provide in a reduced computation time a close approximation of the solution of the functional variational problem. The main principle is to adapt a classical optimization algorithm, the gradient descent method [START_REF] Nocedal | Numerical optimization[END_REF], to the space of probability measures. Based on this principle, we derive here an efficient algorithm for the variational Bayesian framework, based on an optimization in the density probability functions. Our algorithm is based on the exponentiated gradient introduced in [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] for discrete measures. The first contribution of this paper is to set out the definition and the convergence properties of this algorithm in the probability density functions set. In Section 3 we also illustrate its application in the context induced by the variational Bayesian methodology.

The second contribution, exposed in Section 4, consists in the application of the mentioned methodology to linear inverse problems given by the Gaussian white noise model. Concerning the prior distribution, we emphasizeinformation of sparsity. A sparse Bayesian methodology using a Laplace prior was developed by Seeger [38] [40] or Elad [START_REF] Zibulevsky | L1-L2 optimization in signal and image processing[END_REF]. But more generally, a sparse prior information is introduced by heavy-tailed distributions, such as Bernoulli Gaussian law, [START_REF] Goussard | Maximum a posteriori dectection-estimation of Bernoulli-Gaussian processes[END_REF], mixtures of Gaussian [START_REF] Snoussi | Bayesian unsupervised learning for source separation with mixture of gaussians prior[END_REF][START_REF] Ichir | Hidden markov models for wavelet-based blind source separation[END_REF], Cauchy distribution, or α-stable distributions. The Gaussian Scale Mixture class, a generalization of the previous ones, is introduced in [START_REF] Wainwright | Scale Mixtures of Gaussians and the statistics of natural images[END_REF] as a model of wavelet coefficients of natural images. The main advantage of Gaussian Scale Mixtures is that they can moreover be easily written as Gaussian distributions, conditioned by a hidden variable and thus allow the use of Gaussian based tools.

In Section 5, we present simulation results, first on a tomographic example where the a priori information promotes pulses or extremely sparse solutions and secondly on an identification problem in a very large dictionary learning context. Finally, Section 6 concludes the paper.

2. Optimization algorithm in a measures space. In this part, we assume that we stand in the measurable space (R N , B(R N )), where B(R N ) is the σ-field of Borel sets of R N . This section is devoted to the construction of a gradient based algorithm adapted to the optimization problem involved in the variational Bayesian methodology. This problem can be seen as a maximization of a concave functional in a probability density functions set.

Concerning probability density functions there are several possible representations of such objects. The first one is to consider that this set is a subset of the positive cone of L 1 (R N ). As L 1 (R N ) is a Banach space, classical optimization algorithms holds in this case. However, one has to pay a particular attention to the fact that the positivity of the functions together with the fixed total mass imposes additional constraints which can be hardly tractable in some case. Another point of view, see [START_REF] Molchanov | Tangent sets in the space of measures: with applications to variational analysis[END_REF], is to consider this set as a subset of the space of signed Radon measures M(R N ), that is measures that can be written as µ = µ +µ -, endowed with the norm of total variation. Once again this is a Banach space, see [START_REF] Molchanov | Tangent sets in the space of measures: with applications to variational analysis[END_REF]. The classical gradient descent can also be adapted to this framework, as shown in [START_REF] Molchanov | Steepest descent algorithms in a space of measures[END_REF]. However in [START_REF] Molchanov | Steepest descent algorithms in a space of measures[END_REF], measures obtained at each iteration no longer have densities, and this algorithm cannot converge to a solution of our optimization problem. Moreover, this space does not take into account the separability property of measures, which plays a crucial part in our case. Hence in the following we rather consider the Cartesian product of one-dimensional Radon measure spaces, M(R N ) = N i=1 M(R), endowed with the norm:

∀µ ∈ M, µ T V = N i=1 sup Ai∈B(R) Ai dµ + i (x i ) + Ai dµ - i (x i ). (2.1)
Note that when µ is a density measure, i.e. dµ(x) = q(x)dL(x) = q(x)dx, L standing for the Lebesgue measure, its total variation norm coincides with the L 1 norm of its density function q.

Furthermore, a separable probability density function is considered as an element of the closed convex set Ω, defined by

Ω = {µ ∈ M; dµ(x) = N i=1 q i (x i )dx i , where q i ∈ L 1 (R), q i ≥ 0 a.e. and R N dµ(x) = 1}.
(2.2) Note that this set can also be written as the Cartesian product of the sets (Ω i ) i=1,...,N

where

Ω i = {µ i ∈ M(R); dµ i (x) = q i (x i )dx i , where q i ∈ L 1 (R), q i ≥ 0 a.e. and R dµ i (x i ) = 1}.
2.1. Preliminaries. Our purpose in the following is, given a concave functional F : M → R, to define an algorithm which approximates a probability measure µ opt solution of

µ opt = arg max µ∈Ω F (µ). (2.3)
This problem can be seen as a constrained convex optimization problem in the infinite dimensional Banach space ( M, . T V ). In this framework, most results of optimization are based on duality results, see for instance [START_REF] Kellerer | Measure theoretic versions of linear programming[END_REF][START_REF] Chartrand | A Gradient Descent Solution to the Monge-Kantorovich Problem[END_REF]. In the present paper we consider a gradient-like descent algorithm defined directly on the probability measures set, as in [START_REF] Molchanov | Tangent sets in the space of measures: with applications to variational analysis[END_REF].

Let us introduce some notations from [START_REF] Molchanov | Tangent sets in the space of measures: with applications to variational analysis[END_REF]. Let F : M → R. As M is a Banach space, one can compute the Fréchet derivative of F at µ ∈ M as the bounded linear functional dF µ (.) : M → R satisfying

F (µ + ν) -F (µ) -dF µ (ν) = o( ν T V ), when ν → 0.
We also consider the Gateaux derivative of F :

∀ν ∈ M, ∂F µ (ν) = lim t→0 F (µ + tν) -F (µ) t .
In some cases, as in the following one can find a function df : M × R N → R continuous and upper bounded, such that the Gateaux differential of F can be written as

∀ν ∈ M, ∂F µ (ν) = R N df (µ, x)dν(x).
(2.4)

Consider an auxiliary concave functional G : M → R. An important property appearing in the following is that its Fréchet differential is L-Lipschitz on Ω, i.e.

∀(µ 1 , µ 2 ) ∈ Ω 2 , ∀ν ∈ M |dG µ1 (ν) -dG µ2 (ν)| ≤ L ν T V µ 1 -µ 2 T V .
(2.5)

The Lipschitz differential condition of G together with its concavity implies that, see [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] for instance,

∀(µ 1 , µ 2 ) ∈ Ω 2 , 0 ≥ G(µ 1 ) -G(µ 2 ) -dG µ2 (µ 1 -µ 2 ) ≥ -L µ 1 -µ 2 2 T V . (2.6) Furthermore we say that a function F is twice differentiable in the sense of Fréchet at µ ∈ M if dF µ is differentiable. If it exists, d 2 F is a bilinear application from M× M to R.

Statement of the main result.

Let us consider as a first step the hypotheses imposed on the functional F in the rest of this part, given by the following definition.

Hypothesis 1. Let F : M → R be a concave functional. We say that F satisfies Hypothesis 1 if:

(i) F can be written as F = G + H where G is a concave L-Lipschitz Fréchetdifferentiable functional whereas H corresponds to the entropy of the probability measure, that is the Kullback-Leibler divergence from this measure to the Lebesgue measure. (ii) For µ ∈ M, G is twice differentiable in the sense of Fréchet at µ and the first order derivative of F satisfies Equation (2.4). (iii) lim µ →∞ F (µ) = -∞. Remark 1. Hypothesis (i) can be replaced by the more restrictive hypothesis that F is L-Lipschitz Fréchet differentiable.

Our purpose is to construct an iterative algorithm providing a density at each iteration and approximating the solution of (2.3) for a certain class of functionals F . The key principle of our method is given by the Radon-Nikodym Theorem, see [START_REF] Rudin | Real and complex analysis[END_REF] for instance. Let k ≥ 0 be an integer and assume that µ k ∈ M is a probability measure absolutely continuous respectively to the Lebesgue measure. We construct µ k+1 ∈ M as a measure which is absolutely continuous with respect to µ k . In this case, the Radon-Nikodym theorem ensures that this measure should be written as

dµ k+1 (x) = h k (x)dµ k (x), (2.7)
where

h k ∈ L 1 (µ k ) is a positive function. Our aim is to determine a function h k ∈ L 1 (µ k ) which ensures that F (µ k+1 ) ≥ F (µ k
), F being a concave functional satisfying Hypothesis 1. Following the classical iterative scheme given by the gradient descent method, we consider h k as a function of the derivative of F at µ k and, according to our structure

h k (x) = K k (α k ) exp(α k df (µ k , x)), (2.8)
where df is defined by (2.4) whereas α k > 0 is the optimal algorithm step-size at iteration k and K k (α k ) is the normalization constant such that R N dµ k+1 (x) = 1.

We also impose the convention that h k (x) = ∞ when exp(α k df (µ k , x)) is not integrable. One can see that as soon as µ 0 is a measure absolutely continuous with respect to the Lebesgue measure, so is each µ k . This choice of h is motivated by the positive, integrable assumption together with, as mentioned earlier, its coherence with the structure of the gradient descent method. Furthermore the exponential function is the more suitable when dealing with entropy constraint, see [START_REF] Tseng | On the convergence of the exponential multiplier method for convex programming[END_REF] for details. This descent algorithm is defined as the "exponentiated gradient" descent in [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF]. Since [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] it has been widely studied in the context of machine learning even in the Bayesian framework, see [START_REF] Globerson | Exponentiated gradient algorithms for log-linear structured prediction[END_REF] for instance.

The optimization algorithm involved in this paper is the following:

Algorithm 1 Exponentiated Gradient algorithm 1: Initialize( µ 0 ∈ Ω) 2: repeat 3: Compute df (µ k , x) 4: Compute α k = arg max α K k (α) exp(αdf (µ k , .))µ k 5: Compute µ k+1 = K k (α k ) exp(α k df (µ k , .))µ k 6: until Convergence
Let us now state the main convergence result of this paper.

Theorem 2.1. Let F be a concave functional satisfying Hypothesis 1. Let, for every k ≥ 0, α k be the optimal stepsize of Algorithm 1. Then the sequence (µ k ) k≥0 of elements of M given by µ k+1 = K k (α k ) exp(α k df (µ k , .))µ k converges to a maximizer of F on Ω.

2.3. Proof of Theorem 2.1. Before considering the convergence results, let us introduce some differentiability properties of the entropy function in our context. Note that as our purpose is to construct measures with a density, thus we consider in the following that for each k ≥ 0, dµ k = q k dx. The term H in Hypothesis 1 also corresponds in this case to the entropy H of the corresponding q k . Concerning the entropy of density functions, one can notice that it is not differentiable in the sense of Fréchet. However in our case, such as in [START_REF] Frigyik | Functional Bregman Divergence and Bayesian Estimation of Distributions[END_REF], one can find directions h such that the Gateaux differential of H at q k in the direction h can be defined. This is the result stated in the following lemma.

Lemma 2.2. Let us assume that (q k ) k∈N is the sequence of densities corresponding to the measures provided by Algorithm 1. Let k > 0 be fixed and define

W k := {h ∈ L ∞ (q k dL), h > -1 a.s.}.
Then as soon as h ∈ W k , H is differentiable at q k in the direction hq k and

∂H q k (hq k ) = -(ln(q k (x)) + 1)h(x)q k (x)dx.

Proof.

Following the development made in [START_REF] Frigyik | Functional Bregman Divergence and Bayesian Estimation of Distributions[END_REF] for simple functions, let us now prove that when h ∈ W k , the Gateaux derivative of H at q k in the direction hq k can be defined. Hence we consider

H(q k +hq k )-H(q k ) = - R N ln(1+h(x))q k (x)dx- R N ln(q k (x)+h(x)q k (x))h(x)q k (x)dx,
where the existence of each term is ensured by the fact that h ∈ W k .

Let us define

∂H q k (hq k ) = -(ln(q k (x)) + 1)h(x)q k (x)dx, we thus have

H(q k + hq k ) -H(q k ) -∂H q k (hq k ) = - R N ln(1 + h(x))q k (x)dx - R N ln(q k (x))h(x)q k (x)dx -ln(1 + h(x))h(x)q k (x)dx + (ln(q k (x)) + 1)h(x)q k (x)dx) = - R N ln(1 + h(x))(1 + h(x))q k (x)dx + R N h(x)q k (x)dx.
In order to conclude, let us now consider that h is such that h L ∞ → 0. In this case, for almost every x one has h(x) → 0 and ln(1 + h(x)) ∼ h(x). This induces

H(q k +hq k )-H(q k )-∂H q k (hq k ) ∼ - R N h(x)(1+h(x))q k (x)dx+ R N h(x)q k (x)dx = - R N h(x) 2 q k (x)dx,
and as µ k (dx) = q k (x)L(dx) is a probability measure, one has

|H(q k + hq k ) -H(q k ) -∂H q k (h)| h L ∞ ∼ h 2 L 2 (q k dL) h L ∞ (q k dL) ≤ h L 2 (q k dL) ,
which yields the desired result.

Concerning the proof of Theorem 2.1, it involves two main steps. As a first step we prove that the sequence (F (µ k )) k∈N is an increasing sequence. Secondly, we consider the convergence of the sequence (µ k ) k≥0 to a solution of (2.3).

Lemma 2.3. Let F be a functional satisfying Hypothesis 1. Let also (µ k ) k∈N be the sequence provided by Algorithm 1. Then there exist α 0 > 0 such that

∀α ∈ (0, α 0 ) F (µ α ) -F (µ k ) ≥ 0.
(2.9)

Proof.

Let k > 0 be fixed and µ k be given. For every α ≥ 0 we define µ α as the measure defined for all x ∈ R N by dµ

α (x) = K k (α) exp(αdf (µ k , x))dµ k (x) := h α (µ k , x)dµ k (x).
We define furthermore g k (α) := F (µ α ). Thus g k is a function from R + to R and α opt is an optimal step-size if g k (α opt ) = max g k (α), i.e.

α opt = arg max α g k (α).
(2.10)

According to the point (iii) of Hypothesis 1, F (µ) → -∞ when µ → ∞ ensures that we can find an α opt , not necessarily unique, such that ∀α > 0, F (µ α ) ≤ F (µ αopt ).

(2.11)

Let α > 0 be given and consider the decomposition given by point (i) of Hypothesis 1. Thanks to Equation (2.6) one has

G(µ α ) -G(µ k ) ≥ dG µ k (µ α -µ k ) -L µ α -µ k 2 T V .
(2.12)

Furthermore, as µ α = h α (µ k , .)µ k and µ k is a probability measure one can notice that

-L µ α -µ k 2 T V = -L h α (µ k , .) -1 2 L 1 (µ k ) ≥ -L h α (µ k , .) -1 2 L 2 (µ k ) . (2.13)
Furthermore,

H(µ α ) = - R N ln dµ α dL dµ α (x) = H(µ k ) - R N ln(h α (µ k , x))dµ α (x) - R N ln dµ k dL (h α (µ k , x) -1)dµ k (x).
Applying Lemma 2.2 to q k , the density of µ k and h = h α (µ k , .) -1 gives

∂H µ k (µ α -µ k ) = R N (-ln dµ k dL -1)(h α (µ k , x) -1)dµ k (x), (2.14)
This entails

H(µ α )-H(µ k ) ≥ ∂H µ k (µ α -µ k )- R N ln(h α (µ k , x))dµ α (x)- R N (h α (µ k , x)-1)dµ k (x)
(2.15) And thus

H(µ α ) -H(µ k ) -∂H µ k (µ α -µ k ) ≥ - R N ln(h α (µ k , x))dµ α (x), (2.16)
as

R N (h α (µ k , x) -1)dµ k (x) = 0.
Finally, from (2.12), (2.13) and (2.16) one has

F (µ α ) -F (µ k ) ≥ ∂F µ k (µ α -µ k ) -L h α (µ k , .) -1 2 L 2 (µ k ) - R N ln(h α (µ k , x))dµ α (x).
(2.17) Finally, F (µ α ) -F (µ k ) is positive if the right side of Equation (2.17) is positive. This part being quite technical it is included in Appendix 7.1.

Lemma 2.3 ensures that for α > 0 small enough, F (µ α ) ≥ F (µ k ). As we choose µ k+1 = µ αopt , where α opt is defined by (2.10), we obviously have

F (µ k+1 ) ≥ F (µ k ).
Finally the sequence (F (µ k )) k∈N is increasing and upper bounded in R, thus convergent. It thus also satisfies that

F (µ k+1 ) -F (µ k ) → 0.
In order to conclude we have to show that (µ k ) k∈N indeed converges to a maximum of F on Ω. But, for every k ≥ 0, µ k ∈ Ω, which is a bounded set in M and thus in M(R N ). Furthermore, the latter is the dual of C 0 , the space of continuous functions f that tend to zero at infinity, which is a separable Banach space. As a consequence of the Banach-Alaoglu Theorem, see [START_REF] Rudin | Real and complex analysis[END_REF] for instance, there exists µ lim ∈ M and a subsequence (µ kn ) n∈N such that for every continuous function that goes to zero at infinity,

f (x)µ kn (dx) → f (x)µ lim (dx).
i.e. when k → ∞, we have µ kn ⇀ * µ lim ∈ Ω. (µ kn ) n∈N by critical From Lemma 2.3 we know that

F (µ k+1 ) = g k (α opt ) ≥ g k (α), ∀α > 0.
However the analytic form of α opt is not tractable. We thus approximate it by a calculable α subopt , not necessarily smaller than α 0 . In order to determine this α subopt let us determine the differentiability properties of the real valued function g(α) = F (µ α ). From the definition of µ α , it has a density with respect to q k which belongs to W k . As G is assumed to be differentiable, so is F . Furthermore, G is also assumed to be twice differentiable in the sense of Fréchet and from Equation (2.14), so is H in the direction of h α . Hence g k is twice differentiable and the Taylor-Young formula gives, for α small enough,

g k (α) = g k (0) + αg ′ k (0) + α 2 2 g ′′ k (0) + α 2 ε(α) := ϕ k (α) + α 2 ε(α), (2.18)
where ε(α) → 0 when α → 0.

The critical point of ϕ k (α) is α subopt = - g ′ k (0) g ′′ k (0)
, as soon as g ′′ k (0) = 0 exists, which gives in (2.18):

g k (α subopt ) = g k (0) - g ′ k (0) 2 2g ′′ k (0) + α 2 subopt ε(α subopt ), (2.19)
and by construction of µ k+1 ,

F (µ k+1 ) ≥ g k (α subopt ) = g k (0)- g ′ k (0) 2 2g ′′ k (0) +α 2 subopt ε(α subopt ) = F (µ k )- g ′ k (0) 2 2g ′′ k (0) +α 2 subopt ε(α subopt ).
As

F (µ k+1 ) -F (µ k ) → 0, obviously lim k→∞ g ′ k (0) 2 2g ′′ k (0) = 0.
Let us now consider a convergent subsequence (k n ) n∈N and denote by (γ kn ) n∈N the sequence defined ∀n ∈

N by γ kn = - g ′ kn (0) g ′′ kn (0)
, that is the sequence of suboptimal stepsizes. We have then -

g ′ kn (0) 2 g ′′ kn (0) = g ′ kn (0)γ kn → 0.
As df is supposed to be continuous in Equation (2.4), the sequence (g ′ kn (0)) n∈N is convergent and either γ kn → 0 or g ′ kn (0) → 0. Let us assume that γ kn → 0 and that g ′ kn (0) → l = 0. As γ kn → 0 we have, for n large enough,

g kn (γ kn ) -g kn (0) = γ kn g ′ kn (0) + γ 2 kn 2 g ′′ kn (0) + γ 2 kn ε(γ kn ) = 1 2 γ kn g ′ kn (0) + γ 2 kn ε(γ kn ).
Hence,

g kn (γ kn ) -g kn (0) γ kn = 1 2 g ′ kn (0) + γ kn ε(γ kn ), (2.20)
and when n tends to infinity γ kn tends to zero and taking limits in (2.20) one obtains

l = l 2 ,
which is impossible. Hence, g ′ kn (0) → 0 when n → ∞. Let us determine the derivative of the different functions g k . For this purpose, we have to determine the derivative of the function g :

α → h α (µ k , .). As h α (µ k , x) = K k (α)e αdf (µ k ,x) , its derivative is given by ∀x ∈ R N , ∂ g ∂α (α, x) = ∂K k ∂α (α)e αdf (µ k ,x) + df (µ k , x)K k (α)e αdf (µ k ,x) . (2.21)
As α is supposed to be close to zero, one can assume that it is not greater than one and the Lebesgue's Theorem thus allows to invert differentiation and integration in the evaluation of ∂K k ∂α . As

K k (α) = 1 e αdf (µ k ,x) dµ k (x) , (2.22)
we have

∂K k ∂α (α) = - df (µ k , x)e αdf (µ k ,x) dµ k (x) e αdf (µ k ,x) dµ k (x) 2 = -K k (α) df (µ k , x)h α (µ k , x)dµ k (x) = -K k (α) df (µ k , x)dµ α (x).
(2.23)

Combining (2.21) and (2.23) yields ∀x ∈ R N , ∂ g ∂α (α, x)dµ k (x) = dµ α (x) df (µ k , x) - R N df (µ k , y)dµ α (y) .
(2.24) And

g ′ k (α) = ∂F µ α ( ∂ g ∂α µ k ) = R N df (µ α , x) df (µ k , x) - R N df (µ k , y)dµ α (y) dµ α (x),
which leads to

g ′ k (0) = ∂F µ k ( g ′ (0)µ k ) = R N df (µ k , x) df (µ k , x) - R N df (µ k , y)dµ k (y) dµ k (x) = R N df (µ k , x) 2 µ k (dx) - R N df (µ k , y)dµ k (y) 2 (2.25)
Hence, for n large enough,

df (µ kn , .) 2 L 2 (µ kn ) -df (µ kn , .) 2 L 1 (µ kn ) → 0,
and df (µ kn , .) converges to a constant λ, independent of x ∈ R N . Let ν be any element of Ω, then ∂F µ lim (νµ lim ) = 0, and from concavity of F we have

F (ν) ≤ F (µ lim ) ∀ν ∈ Ω.
Which means that F (µ lim ) is a maximum of F over Ω.

In the present part we have presented the convergence properties of our algorithm, which is well adapted to the space of probability measures. Our interest in the following is also in its application in the context of variational Bayesian methodology. For the sake of completeness, let us recall this methodology introduced in [START_REF] Mackay | Ensemble learning and evidence maximization[END_REF].

Application to variational Bayesian methodology.

3.1. Classical variational Bayesian algorithm. For the sake of simplicity we consider in the following density functions q instead of the corresponding measures µ. We denote by y ∈ R M the M dimensional vector containing the data information whereas w ∈ R N represents the vector to be estimated, which is considered as a realization of a random vector W . We also denote by p the prior probability density function (p.d.f.) of W . The Bayes rule entails that this prior distribution is closely related to the posterior one, p(w|y), up to a normalization constant. Even in simple cases this posterior may not be separable. Hence, in the variational Bayesian framework, we approximate it by a separable probability density q(w) = i q i (w i ).

(3.1)

Taking separable laws obviously simplifies the problem even if it introduces some approximation errors.

Therefore, the purpose of variational Bayesian methodology is to determine a separable probability density function q close to the true posterior in the sense defined by the Kullback-Leibler divergence see [START_REF] Smidl | The variational bayes method in signal processing[END_REF] for instance.

The optimal approximating density q is then given by ∀i ∈ {1, . . . , N },

q i (w i ) = K i k exp < ln p(y, w) > j =i qj , (3.2)
where K i k is the normalization constant and < ln p(y, w) > j =i qj = R N -1 ln p(y, w) j =i q j (w j )dw j is the mean of ln p(y, w) under the probability j =i q j . Although this solution is obtained analytically, Equation (3.2) clearly does not have an explicit form. In order to have implementable methods, several approaches can be considered. The first is to impose conjugate prior distributions which ensure that the posterior distribution belongs to a given family and to reduce the optimization problem to an approximation of its parameters. A different approach, introduced in [START_REF] Seeger | Large scale bayesian inference and experimental design for sparse linear models[END_REF] is to consider approximating law in a given family, namely the Gaussian one, which does not necessarily contains the posterior distribution. These approaches also reduce the functional optimization problem to a parametric one. However even in these cases the intricate form of (3.2) imposes a numerical approximation of the solution. A natural method in this case is to consider alternate descent algorithms. However this method is time consuming, each update step needing N iterations, one for each component q i . This drawback can be diminished by relaxing the separability assumption and considering instead a descent by groups of coordinates. In this case the main drawback is given by the covariance matrix of each group which needs to be stored and inverted. Several techniques, such as the numerical preprocessing of this matrix in [START_REF]Variational bayesian super resolution[END_REF] have thus been considered to overcome this drawback. In [START_REF] Seeger | Large scale bayesian inference and experimental design for sparse linear models[END_REF] an alternative way is developed, based on concave-convex programming. But for very large dimensional problems these methods remain inefficient in general, as they need several intricate approximations. In the present work we consider instead the application of the algorithm defined in Section 5.1.1 to the functional optimization problem induced by the Bayesian variational problem.

Variational Bayesian Exponentiated Gradient Algorithm.

In this section we define an iterative method which allows to compute efficiently at each iteration each q i independently of the other ones.

A first step is to rewrite the minimization problem as a convex optimization problem independent of the posterior distribution to be approximated. Instead of minimizing the Kullback-Leibler divergence, we thus remark, as in [START_REF] Choudrey | Variational methods for bayesian independent component analysis[END_REF], that ln p(y) = ln p(y, w) p(w|y) ,

where w is the vector of hidden variables and parameters.

As the log-likelihood ln p(y) in (3.3) does not depends on w one can write

ln p(y) = F(q) + KL[q||p(.|y))].
In this case,

F(q) = R N q(w) ln p(y, w) q(w) dw, (3.4)
is the negative free energy. Thus, minimizing the Kullback-Leibler divergence is obviously equivalent to maximizing this negative free entropy.

Therefore, in the following we consider the problem of maximizing

F(q) = R N ln p(y, w)q(w)dw - R N ln(q(w))q(w)dw = ln p(y, .) q + H(q), (3.5)
where

H(q) = - R N ln(q(w))q(w)dw,
is the entropy of q. The main advantage of this approach is that the objective functional does not depend on the true posterior anymore but only on the joint distribution p(y, .), which is more easily tractable.

One can also notice that the problem of finding

q opt = arg max q separable p.d.f. F(q) (3.6)
is equivalent to the problem of finding

µ opt = arg max µ∈Ω F (µ). (3.7)
Where the functional F is defined F (µ) = F(q), as soon as µ ∈ Ω is such that q is the density of µ. Let us also define in the following ∂f by ∀x ∈ R N , ∂f (q, x) = df (qL, x).

A classical method in this context is to consider each density q as a L 1 (R N ) function and to apply classical optimization algorithms. In the present framework, taking the non-negativity and the total mass assumptions into account, the algorithm involved is given by the projected gradient method which gives:

∀w ∈ R N q k+1 (w) = P Θ (q k (w) + ρ k ∂f (q k , w)), (3.8)
where P Θ is the projector operator on the subspace Θ = {q ∈ L 1 (R N ); q(w) ≥ 0 and q L 1 = 1} and ∂f denotes the Fréchet differential of F at q k . However, this algorithm requires that ∂f (q k , w) ∈ L 1 (R N ) which is not the case in general.

Therefore in the following, we rather apply Theorem 2.1 and the algorithm introduced in Section 2 to the Variational Bayesian framework of Section 3.1. One can easily see that the function F can be written as an entropy term added to a differentiable function G = ln p(y, .) q . As ln p(y, .) does not depend on q, this function is clearly Lipschitz differentiable respectively to q and twice differentiable in the sense of Fréchet. Hence F satisfies Hypothesis 1.

We consider F(q) = ln p(y, .) q + H(q). In this case, the Gateaux differential of F (µ) = F(q) at µ ∈ Ω separable is given by dF

µ (ν) = i R N df (µ i , x i )ν i (dx) with df (µ i , x i ) = ∂f (q i , x i ) such that ∀i ∈ {1, . . . , N }, ∀w ∈ R N , ∂f (q i , w i ) = ln p(y, w) j =i qj -ln q i (w i ) -1.
Let k ≥ 0 be given and q k be constructed. Following the scheme defined by Algorithm 1 and Equation (2.8), at the following iteration we consider q α given, for α > 0, by

∀w ∈ R N , q α k (w) = K k (α k )q k exp α k ∂f (q k , w) (3.9) = Kk (α k )q k (w)   i exp < ln p(y, w) > j =i q k j (wj ) q k i (w i )   α k = Kk (α k )q k (w) i q r i (w i ) q k i (w i ) α k (3.10)
where Kk (α k ) is the normalization constant and q r is an intermediate function defined by ∀i ∈ {1, . . . , N }, q r i (w i ) = exp < ln p(y, w) >

j =i q k j (wj )
The main challenge is to determine the value of α k > 0. This optimal value α opt should satisfy g ′ k (α opt ) = 0. However, this quantity is hardly tractable in practice. Therefore, we consider instead the suboptimal value given by

α subopt = - g ′ k (0) g ′′ k (0) , (3.11)
when g ′′ k (0) = 0. This leads to the main algorithm of this paper.

Algorithm 2 Variational Bayesian Exponentiated Gradient Like Algorithm

1: Initialize( q 0 ∈ Θ) 2: repeat 3: function Iteration( Compute q k+1 = K k q k exp α k ∂f (q k , w) ) 4:
Compute q r i (w i ) = exp < ln p(y, w) > j =i q k j (wj ) for every i = 1, . . . , N

5:

Compute q α (w) = Kk (α)q k (w) q r (w) q k (w) α .

6:

Compute α subopt = -

g ′ k (0) g ′′ k (0)
7:

Take q k+1 = q α subopt .

where y ∈ R M is the vector of observations given as a linear function of the unknown vector x ∈ R N to be estimated. Here, b ∈ R M is the noise vector whereas H is a matrix in M N ×M . We also suppose that x is a realization of a random vector X.

In the following we stand in a white noise model which induces that the noise is supposed to be an iid Gaussian vector N (0, σ 2 b I). The corresponding likelihood is

p(y|x) = (2πσ 2 b ) -M/2 exp - y -Hx 2 2σ 2 b . (4.2)
Concerning the prior distribution we choose to take sparsity into account by considering X distributed following a separable heavy tailed distribution. The most general case is given by Gaussian Vector Scale Mixture (GVSM) defined in [START_REF] Jana | Optimality of KLT for high-rate transform coding of gaussian vectorscale mixtures : Application to reconstruction, estimation, and classification[END_REF]. In this case, for i = 1, . . . , N , we suppose that X i ∼ U i / √ Z i where U ∼ N (0, σ 2 s I), Z = Z i is a positive random vector of independent positive coordinates and U and Z are independents. As a consequence the density of X is given in an integral form as

∀i ∈ {1, . . . , N }, p(x i ) = R √ z i (2π) 1/2 σ s e - z i x 2 i 2σ 2 s φ zi (z i )dz i .
Note that in the definition, for the sake of simplicity, we consider Z as a precision parameter and not a dilatation one. Gaussian Vector Scale Mixture forms a large class of nongaussian random variables recently developed as a model of wavelet coefficients of natural images, see [START_REF] Wainwright | Scale Mixtures of Gaussians and the statistics of natural images[END_REF]. The main interest of this model is, by solving an extended problem due to the presence of a hidden random vector Z, to allow the use of Bayesian hierarchic approaches. The Gaussian Scale Mixture family offers a large class of random variables including Gaussian mixture, when Z ∼ ZI a discrete random vector or Student laws if the Z i correspond to Gamma random variables. With different hypothesis on the distribution of Z one can also define Generalized Gaussian distributions or α-stable ones, see [START_REF] Wainwright | Scale Mixtures of Gaussians and the statistics of natural images[END_REF]. Indeed GSM offers a simple representation of a large class of nongaussian probability distributions, which justify the increasing interest in this model.

In our context, we choose to consider Z as a independent Gamma random vector, i.e. for i = 1, . . . , N , we have

Z i ∼ G(ã i , bi ) and ∀i ∈ {1, . . . , N }, p(x i ) = bãi i Γ(ã i ) R √ z i (2π) 1/2 σ s e - z i x 2 i 2σ 2 s z ãi-1 i e -zi bi dz i . (4.3)
For ãi = bi = ν 2 the p.d.f. of X corresponds to a Student-t distribution, as in the model used in [START_REF] Chantas | Variational Bayesian image restoration based on a product of t-distributions image prior[END_REF]. This model of Z ensures that X satisfies the conjugate priors condition.

One can easily check that when the prior information is given by (4.3), Equation (4.2) gives the following posterior distribution

p(x, z|y) ∝ σ -M b exp - y -Hx 2 2σ 2 b N i=1 √ z i σ s exp - z i x 2 i 2σ 2 s bãi i z ãi-1 i e -zi bi Γ(ã i ) . (4.4)
Considering that we do not know the constants involved and that the mixing matrix H is high dimensional, this posterior distribution cannot be evaluated directly.

4.2. Numerical implementation. The aim of variational Bayesian methodology and therefore of our method in the context established in Part 4 is the approximation of the posterior p.d.f given by (4.4) by a separable one q(x, z) = i q i (x i ) j qj (z j ). As we have chosen conjugate priors for X and Z, the optimum approximating distribution of X is known to belong to a Gaussian family, whereas the p.d.f. of Z belongs to a Gamma one.

q k ∼ i N (m k (i), σ 2 k (i)) qk ∼ j G(a k (j), b k (j))
Hence at the initialization stage, we consider

q 0 = N (m 0 , Diag(σ 2 0 )) q0 = j G(a 0 (j), b 0 (j))
where Diag(v) is a diagonal matrix with v on its diagonal, and σ 2 0 ∈ R N is the vector of initial variances.

Our minimization problem can be analyzed following the alternate iterative scheme: qk+1 = arg max q F q k q q k+1 = arg max q F q qk+1 4.2.1. Approximation of q. One can see in Equation (4.4) that the conditional posterior p(z|x, y) is fully separable. In this case the classical Bayesian variational approach is efficient enough to be implemented directly. Hence all the qk+1 i can be computed simultaneously, knowing only q k . Thanks to the classical variational Bayesian approach, [START_REF] Mackay | Ensemble learning and evidence maximization[END_REF], described in Section 3, we deduce qk+1 thanks to Equation (3.2) and Equation (4.4), for every i = 1, . . . , N

qk+1 i (z i ) ∝ exp < ln p(y, x, z) > j =i qk j (zj )q k (x) (4.5) ∝ exp   (ã i - 1 2 ) ln(z i ) -( x 2 i z i 2σ 2 s + z i bi ) l q k l (x l ) j =i qk j (z j )dxdz   ∝ exp (ã i - 1 2 ) ln(z i ) -z i bi - x 2 i z i 2σ 2 s q k i (x i )dx i ∝ exp (ã i - 1 2 ) ln(z i ) -z i bi - (σ 2 k (i) + m 2 k (i))z i 2σ 2 s ∝ z ãi-1 2 i exp -z i bi + (σ 2 k (i) + m 2 k (i)) 2σ 2 s (4.6)
This entails that qk+1 i corresponds to a Gamma p.d.f. of parameters: ∀i ∈ {1, . . . , N },

a k+1 (i) = ãi + 1 2 , (4.7) b k+1 (i) = m 2 k (i) + σ 2 k (i) 2σ 2 s + bi . (4.8)
4.2.2. Approximation of q by Algorithm 2. Let us assume that at the initialization stage, q 0 is a Gaussian p.d.f. with mean m 0 and covariance matrix Diag(σ 2 0 ). At each iteration k +1 we determine the approximation of q k+1 thanks to our method. Let us define an auxiliary function q r by q r i (x i ) = exp ln p(y, x, z)

j =i q k j (xj )q k+1 (z) , thus ∀i ∈ {1, . . . , N }, q r i (x i ) = exp < ln p(y, x, z) > j =i q k j (xj )q k+1 (z) (4.9) ∝ exp   - y -Hx 2 2σ 2 b + x 2 i z i 2σ 2 s j =i q k (x j )q k+1 (z)dxdz   ∝ exp   - x T H T Hx -2x T H T y 2σ 2 b + x 2 i z i 2σ 2 s j =i q k j (x j )q k+1 (z i )dxdz i   ∝ exp - 1 2σ 2 b x 2 i diag(H T H) i -2x i (H T y) i + 2x i (H T Hm k ) i -2x i diag(H T H)) i m k (i) + x 2 i a k+1 (i) 2σ 2 s b k+1 (i) (4.10)
where diag(A) is the vector composed by the diagonal entries of A. Note that q r corresponds, up to the normalization term, to the density of a Gaussian distribution with mean m r and variance σ 2 r , where, for every i = 1, . . . , N ,

σ 2 r (i) = diag(H T H) i σ 2 b + a k+1 (i) b k+1 (i)σ 2 s -1 (4.11)
and

m r (i) = σ 2 r (i) × H T y -(H T H -diag(H T H))m k σ 2 b i (4.12)
Therefore, by Equation (3.10), we have for every i = 1, . . . , N ,

q α i (x i ) = K k (α)q k i (x i ) q r i (x i ) q k i (x i ) α = σ 2 k (i) σ 2 r (i) α K k (α) exp - (x i -m k (i)) 2 2σ 2 k (i) exp -α x 2 i (σ 2 k (i) -σ 2 r (i)) 2σ 2 r (i)σ 2 k (i) × exp -α -2x i (m r (i)σ 2 k (i) -m k (i)σ 2 r (i)) + m r (i) 2 σ 2 k (i) -m k (i) 2 σ 2 r (i) 2σ 2 r (i)σ 2 k (i) = σ 2 k (i) σ 2 r (i) α K k (α) exp - 1 2 x 2 i σ 2 r (i) + α(σ 2 k (i) -σ 2 r (i)) σ 2 r (i)σ 2 k (i) × exp - 1 2 -2x i m k (i)σ 2 r (i) + α(m r (i)σ 2 k (i) -m k (i)σ 2 r (i)) σ 2 r (i)σ 2 k (i) + t(α)
where q α is defined in Section 5.1.1, and t(α) = α

mr(i) 2 σ 2 k (i)-m k (i) 2 σ 2 r (i) 2σ 2 r (i)σ 2 k (i)
is a constant. Finally, q α i still corresponds to a Gaussian p.d.f. with parameters m α and Diag(σ 2 α ) satisfying:

σ 2 α (i) = σ 2 r (i)σ 2 k (i) σ 2 r (i) + α(σ 2 k (i) -σ 2 r (i)) (4.13) m α (i) = m k (i)σ 2 r (i) + α(m r (i)σ 2 k (i) -m k (i)σ 2 r (i)) σ 2 r (i) + α(σ 2 k (i) -σ 2 r (i)) . (4.14)
In order to construct q k+1 we choose in the previous equation α = α subopt defined in Equation (3.11).

Finally, we obtain the following algorithm.

Algorithm 3 Supervised Sparse Reconstruction algorithm (SSR)

1: Initialize(q 0 , q0 ) 2: repeat function Estimate q k+1 (q k+1 ) 8: compute q r ← (m r , σ 2 r ) by Equation (4.12) and Equation (4.11) compute α subopt and Equation (??).

9: compute q α ← (m α , σ 2
11:

compute q k+1 = q α subopt .

12:

end function 13: until Convergence 4.3. Unsupervised algorithm. The algorithm described in the previous part is not a fully Bayesian one as it still depends on some hyperparameters, namely the variances induced by the model (4.1) and (4.3). We see in the following how this method can be extended to an unsupervised one by estimating these parameters. The parameters of the underlying Gamma random variable are not estimated in the following as they define the sharpness of the prior distribution. We thus only estimate the variance parameter of this prior together with the trade off between the prior and the noise.

In order to simplify the different expressions, we introduce in the following the notations γ b = 1/σ 2 b and γ s = 1/σ 2 s . Hence, γ b and γ s are the precision parameters of the noise and of the prior distribution. From now on they are also assumed to be random variable with Gamma prior of parameters (ã b , bb ) resp. (ã s , bs ). As we do not have any information on these precision parameters γ b and γ s , this prior is a Jeffrey's prior obtained by fixing ãb = 0, bb = 0 resp. ãs = 0, bs = 0.

With these assumptions, the posterior distribution from Equation (4.4) can be written as

p(x, z, γ b , γ s |y) ∝ γ M 2 b exp - γ b y -Hx 2 2 γ N 2 s i √ z i exp - γ s z i x 2 i 2 bãi i z ãi-1 i e -zi bi Γ(ã i ) × bã b b γ ãb -1 b e -γ b bb Γ(ã b ) bãs s γ ãs-1 s e -γs bs Γ(ã s ) . (4.15)
As in the previous section, the conditional posterior p(z, γ b , γ s |x, y) is separable and can be approximated thanks to the classical variational Bayesian approach. Once again only the distribution of X needs the use of Algorithm 2. Here the alternate optimization scheme to carry out is:

qk+1 = arg max q F q k qq k b q k s q k+1 = arg max q F q qk+1 q k b q k s q k+1 b = arg max q b F q k+1 qk+1 q b q k s q k+1 s = arg max qs F q k+1 qk+1 q k+1 b q s 4.3.1.
Optimization of the approximate p.d.f. q b . Concerning the random vectors Z and X, the updating process follows the same scheme as the supervised case, see Section 4.2, and is not recalled here. The main differences reside in the update of the parameter distributions.

As the distributions of γ b and γ s are supposed to be Gamma, which is coherent with the conjugate priors hypothesis, at each iteration we just adapt the parameters. Hence we initialize our algorithm by considering that

q 0 b ∼ G(a 0 b , b 0 b )
At iteration k + 1 we consider the maximum of the free energy from Equation (3.2) which gives

q k+1 b (γ b ) ∝ exp < ln p(x, y, z, γ b , γ s ) > j q k+1 j (xj ) j qk+1 j (zj )q k s (γs) ∝ exp < ( M 2 + ãb -1) ln(γ b ) -γ b ( y -Hx 2 2 + bb ) > j q k+1 j (xj ) ∝ exp ( M 2 + ãb -1) ln(γ b ) -γ b 1 2 y -Hm k+1 2 + 1 2 N i=1 diag H t H i σ 2 k+1 (i) + bb So q k+1
b is a Gamma p.d.f. of parameters:

a k+1 b = M 2 + ãb (4.16) b k+1 b = 1 2 y -Hm k+1 2 + 1 2 N i=1 diag H t H i σ 2 k+1 (i) + bb (4.17)
4.3.2. Optimization of the approximate p.d.f. q s . As for γ b , the approximation of q s is performed in the family of Gamma distributions. Hence at the initialization step we assume that

q 0 s ∼ G(a 0 s , b 0 s )
and at iteration k + 1, thanks again to Equation (3.2), we obtain

q k+1 s (γ s ) ∝ exp < ln p(x, y, z, γ b , γ s ) > j q k+1 j (xj ) j qk+1 j (zj )q k+1 b (γ b ) ∝ exp < ( N 2 + ãs -1) ln(γ s ) -γ s ( 1 2 N i=1 z i x 2 i + bs ) > j q k+1 j (xj )q k+1 j (zj ) ∝ exp < ( N 2 + ãs -1) ln(γ s ) -γ s ( 1 2 N i=1 a k+1 (i) b k+1 (i) (m 2 k+1 (i) + σ 2 k+1 (i)) + bs ) So q k+1
s is a Gamma p.d.f. and its parameters are deduced by identification.

a k+1 s = N 2 + ãs (4.18) b k+1 s = 1 2 N i=1 a k+1 (i) b k+1 (i) (m 2 k+1 (i) + σ 2 k+1 (i)) + bs (4.19)
Finally, the algorithm can be summed up as follows.

Algorithm 4 UnSupervised Sparse Reconstruction algorithm (USSR)

1: Initialize(q 0 , q0 , q 0 b , q 0 s ) 2: repeat 3:

function Estimate qk+1 (q k , q k b , q k s )

4:

update a k+1 using Equation (4.7)

5:

update b k+1 using Equation (4.8)

6:
end function

7:
function Estimate q k+1 (q k+1 , q k b , q k s )

8:

compute q r (x) ← (m r , σ 2 r ) using Equation (4.12) and Equation (4.11) compute q k+1 = q α subopt .

9: compute q α (x) ← (m α , σ 2
12:

end function

13:
function Estimate q k+1 b (q k+1 , q k+1 ) end function 21: until Convergence 5. Simulations. This section is devoted to numerical validations of the method proposed in this paper. For the sake of completeness we will treat two inverses problems. The first one is given by a tomographic example. The goal is to enhance the accuracy and the effectiveness of our approach, in comparison with classical ones, such as classical Variational Bayesian methods or Monte Carlos Markov Chain (MCMC) methods. From the limitations of these concurrent approaches, we choose to consider only a small dimensional inverse problem (4096 unknowns), and thus to invert the Radon transform of a small sparse image (64 × 64 pixels).

The second experimental result is devoted to a relatively large inverse problem (≈ 300000 unknowns). In this case, the problem is to identify different components in a dictionary learning process. This learning is performed in a very noisy environment, where the signal-to-noise ratio can have negative values. This signal processing problem can appear for instance in an astrophysical context (detection of gravitational waves [START_REF] Pai | Best network chirplet chain: Near-optimal coherent detection of unmodeled gravitational wave chirps with a network of detectors[END_REF]) or in radar imaging [START_REF] Wang | Inverse synthetic aperture radar imaging of maneuvering targets based on chirplet decomposition[END_REF][START_REF] Amiri | Chirplet based denoising of reflected radar signals[END_REF]. This second example illustrates the performances of the method for large dimensional problems.

In both cases the sparsity information is introduced by an iid Student-t prior. This prior is a particular case of GVSM. In the following we thus fix ãi = ν 2 and bi = ν 2 in Equation (4.4). 5.1. Tomographic example. For the sake of completeness, a short description of the concurrent approaches is given, enhancing the main differences between them. In a second part, we describe the phantom together with the simulation parameters.

Algorithms descriptions.

Filtered Back Projection (FBP). Filtered Back Projection is the classical approach to invert the Radon transform [START_REF] Natterer | Algorithms in tomography, The State of the Art in Numerical Analysis[END_REF][START_REF] Kak | Principles of computerized tomographic imaging[END_REF]. This algorithm is obtained by sampling the continuous inversion formula. Each line of the sinogram (see Fig. 5.1) is filtered with a ramp filter. The filtered data are backprojected. The discrete version of the backprojection operator is given by H t .

Monte Carlos Markov Chain. The MCMC method contains a large class of Bayesian algorithms [START_REF] Robert | Monte-Carlo statistical methods[END_REF]. In the following we consider the Gibbs algorithm for its efficiency when the size of the problem increases. This algorithm is used for a wide class of applications [START_REF] Makni | Joint detection-estimation of brain activity in functional MRI: a multichannel deconvolution solution[END_REF][START_REF] Giovannelli | Unsupervised bayesian convex deconvolution based on a field with an explicit partition function[END_REF][START_REF] Bidon | Knowledge-aided stap in heterogeneous clutter using a hierarchical bayesian algorithm[END_REF]. The principle is to obtain samples of the posterior law given by Equation (4.4) by an alternate sampling with conditional laws. The algorithm is as follows:

(i) z k is sampled with p(z|y, x k-1 ) (ii) x k is sampled with p(x|y, z k ) (iii) go to i) until convergence of the Markov chain. As the conditional law p(z|y, x k-1 ) is a separable Gamma distribution, the computation of the sample z k is easy. Furthermore p(x|y, z k ) is a correlated Gaussian distribution with a covariance matrix

R k = M t k M k = [ 1 σ 2 b H t H + 1 σ 2 s Diag(z k )] -1 and a mean m k = 1 σ 2 b R k H t y.
The sampling under this correlated distribution is performed by sampling a vector of centered iid Gaussian random variables with variance 1. This vector is then multiplied by the correlation matrix M k and added to m k .

Remark 2. At each sampling iteration the covariance matrix of size N × N has to be inverted.

Classical Bayesian Variational approach. This method was already described in Section 3. In the alternate descent algorithm, one can choose the degree of separability of the approximate distribution. In the following we consider two cases. In the first case, the so called VBBloc, we consider, as in [START_REF]Variational bayesian super resolution[END_REF] that the separation of the approximating law is only between x and z. This leads to a consideration of the approximating distribution as: q(x, z) = q(x)q(z)

and q ∼ N (m, R) q ∼ G(a, b)
Thus, with Equation (3.2), we obtain ∀i ∈ {1, . . . , N } the following updating equations, at iteration k + 1:

a k+1 (i) = ν 2 + 1 2 , b k+1 (i) = ν 2 + m 2 k (i) + diag(R k )(i) 2σ 2 s R k+1 = 1 σ 2 s Diag(a/b k+1 ) + 1 σ 2 b H t H -1 m k+1 = 1 σ 2 b R k+1 H t y.
Remark 3. At each step, the updating of the covariance matrix requires the inversion of an N × N matrix, but the convergence rate is better than for the MCMC approach.

To overcome the limit given by a matrix inversion in the classical variational Bayesian framework, we can construct an approximate distribution separable on X. Hence, we estimate a vector of variances instead of the matrix of covariance. This approach is called VBComp in the following.

q(x, z) = Π i q i (x i )q(z)
In this case Equation (3.2) gives the following updating equations, ∀i ∈ {1, . . . , N }:

a(i) k+1 = ν 2 + 1 2 , b k+1 (i) = ν 2 + m 2 k (i) + σ 2 k (i) 2σ 2 s
And, for every i ∈ {1, . . . , N }:

σ 2 k+1 (i) = 1 σ 2 s a(i)/b k+1 (i) + 1 σ 2 b (H t H) (i,i) -1 m k+1 (i) = σ 2 k+1 (i) σ 2 b (H t y(i) -(d(i) -(H t H) (i,i) m k (i))) d = H t Hm k Remark 4.
For each pixel x i , the corresponding value of d = H t Hm k must be determined.

Simulation configuration.

The test image is given by a sparse phantom, composed of 7 peaks on a grid 64 × 64 (see Table 5 Table 5.1 Peaks definition in the phantom been simulated in a parallel beam geometry. These projections are collected from 32 angles θ, uniformly spaced over [0, 180[. Each projection is composed of 95 detector cells. We add a white Gaussian noise (iid) with a standard deviation equal to 0.3 (see Fig. 5.1). Data have thus a relatively bad signal to noise ratio and the number of unknowns is larger than the number of data, which leads to an ill-posed inverse problem. 5.1.3. Results and discussion. In this section, we expose the inversion of this tomographic problem with the approaches described earlier.

All the iterative approaches are initialized with a zero mean and a variance equal to one, and the hyperparameters σ 2 b , σ 2 s and ν are respectively fixed at 1, 0.05 and 0.1. The original image and its different reconstructions given by Filtered BackProjection and posterior mean estimators, are summed up on Fig. 5.2. A comparison of Fig. 5.2 (b) with 5.2 (c), 5.2 (d) and 5.2 (e) clearly shows that the analytical inversion of the Radon transform performed by Filtered Back Projection (FBP) is less robust to noise than Bayesian approaches. Asymptotically, in Bayesian cases, theoretical results are favorable to the MCMC approach as they do not need any approximation. In practice, the number of samples is too small to fit with the asymptotic results of MCMC methods, which explains the bad reconstruction observed in Fig. 5.2(c). Finally, the Supervised Sparse Reconstruction (SSR) (see Fig. 5.2(f)) has the same reconstruction quality as the classical variational Bayesian approaches (see VBBloc Fig. 5.2(d) and VBComp Fig. 5.2(e)). However, when we compare the execution time (see Table 5.2), we see that our approach is 10 time faster than the VBBloc approach, 40 time faster than the VBComp approach and 370 faster than the MCMC approach for this small inverse problem. Moreover, this ratio increases with the size of the problem as both MCMC and VBBloc need the inversion of a covariance matrix at each iteration, which is not the case for our algorithm. 5.1.4. Hyperparameters estimation. As seen in Section 4.2, our approach is defined in a fully Bayesian framework. We thus estimate the values of hyperparameters The results are summed up in Fig. 5.3. We observe that the hyperparameters estimation intensifies the sparsity of the reconstructed image together with the snr, as it goes from 6.00 db in the previous case to 10.06 db. Estimating the true hyperparameters is in this case more relevant than arbitrarily chosen parameters. We observe on Fig. 5.3 (c) that the background is equal to zero even if some additional peaks appear in the reconstructed image.

Finally we see in Fig. 5.3, (d) and (e), the behavior of the estimation of the hyperparameters respective to the number of iterations. This plot is in a logarithm scale due to the dynamic of the fluctuations. We observe that for σ 2 b this estimation converges to the true value (dashed line). For σ 2 s we do not known the true value, but one can notice that this estimation also converges.

Component identification (denoising).

The main advantage of our approach is that it can address larger problems (≈ 300000 unknowns in the present case). In the present part, our purpose is to treat an identification problem in a redundant dictionary decomposition such as in [START_REF] Wipf | Sparse Bayesian learning for basis selection[END_REF][START_REF] Buchgraber | A sliding-window online fast variational sparse Bayesian learning algorithm[END_REF][START_REF] Shutin | Fast adaptive variational sparse Bayesian learning with automatic relevance determination[END_REF]. More precisely, we identify chirps components in a linear mixture. Unfortunately, this mixture is corrupted by noise and spurious signal. To identify each chirp in the mixture and to remove the effect of the spurious signal, we develop a dictionary learning approach. For the sake of simplicity, the construction of the dictionary is included in Appendix 7.2.

The simulated data are sampled at a frequency F e = 44khz, and they are composed of N = 2 16 points, thus the duration of measurement T mes is equal to 1.5 second. The dictionary components are stored in a matrix H and we denote by x ∈ R N the vector of weights to be estimated. Hence, the measurements can be modeled as

y = Hx + b, (5.1)
where b ∼ N (0, σ 2 b I) is a Gaussian white noise. Note that, as described in Appendix 7.2, H represents here an overcomplete dictionary.

We also consider different simulation contexts, with different noise levels and a growing number of components. In each case the simulated data are composed of N = 2 16 points.

As a first step, the data are composed of two cosine terms together with height chirp functions whose parameters are given in Table 5.3. Simulated data have a relatively bad snr, at 5.68 db. These data are given in Fig. 5.4. This inverse problem is then solved thanks to our unsupervised approach. Our algorithm was launched with the shape parameter of the Student-t equals to ν = 0.01 in order to introduce a very sparse prior. The initialization parameters are:

• The mean of q 0 , m 0 = 0,

• the variance of q 0 , σ 2 0 = 1, • the mean of q 0 b , is equal to 10 -5 , • the mean of q 0 s , is equal to 10 5 . After 400 iterations (316 s), the algorithm (USSR) converges to a stable solution. We consider here that the estimation x is obtained by the mean of the approximate distribution q. Our method find the right position of the coefficients in the signal. Concerning the amplitude, the relative error obtained is between 0.4% and 2%, thus very close to the right value. Estimation ŝ of s is performed thanks to the estimation of the coefficients x and Equation (7.6). The snr of ŝ is equal to 22.6 db.

As a second step, we show the performances of the algorithm USSR when the signal s is hidden by the noise, the snr of data being equal to -5 db (see Fig. 5.5 (a)). We generate a signal composed of four chirps of parameters summed up in Table 5.4. The estimation of coefficients is performed using our USSR algorithm with the same initialization as in the previous case. After 400 iterations we obtain the coefficients illustrated in Fig. 5.5, (c) and (d). The reconstructed coefficients are represented by a line and the true values of the coefficients are marked by a cross. We observe that all coefficients are in the right place and that the peak amplitudes are systematically underestimated, but the estimated values are relatively close to the true ones. Fig. 5.5 (b) points out the estimator ŝ. We see that the shape of the signal is close to the true one, and that when the signal is missing, between 1 and 1,45 s, the level of the residual noise is relatively low. In Table 5.5 we compare the reconstructed signal ŝ with the reconstruction obtained with three classical denoising approaches, namely the Wiener filter, the soft wavelet shrinkage, with the four vanishing moments symmlet, and the hard shrinkage in the STFT domain. In these three methods, we furthermore have to tune a parameter which is the correlation power for the Wiener filter and the threshold for the soft wavelet and the hard STFT shrinkage.We therefore choose the value of this hyperparameter which minimizes the Signal-to-Noise Ratio (snr), which depends on the signal. Our approach hugely increases the snr (20 db) as the noise level is divided by 100 whereas the classical methods reduce the noise only by a factor of 4, 5 and 26 respectively. As expected, Wiener and wavelet approaches are not adapted to the considered problem and give insufficient results. Moreover the USSR approach increases the SNR of 6 db respectively to the redundant decomposition given by STFT shrinkage, even if the STFT decomposition seems well adapted in this context. This result comes from the sparsity of the reconstructed signal, which is larger in the case of the USSR reconstruction. Finally, we study the robustness of the USSR algorithm with respect to the Signal to Noise Ratio. This simulation allows a better understanding of the reconstruction properties of our method. Hence we consider simulated data with 6 different snr (-5, -2, 1, 2, 5, 10) and with low (4) and high ( 16) numbers of components. For each snr and each number of components, we consider 30 datasets, the components of the signal being randomly chosen in the dictionary. We reconstruct these 360 sets of data by the USSR algorithm taking the same configuration and the same initialization as in previous cases.

The results are summed up in At first sight we see that there are no false positives with our approach. Indeed, as the approach is unsupervised with a sparsity prior, the coefficients with low energy are considered as noise. Moreover we can reconstruct 16 components without loss when the snr is greater or equal to 5 db and resp. 4 components when the snr is greater to 1 db. Fig. 5.6 shows that the reconstruction is more difficult when the number of components increases. Fig. 5.6 (c) is obtained by calculating the snr of the reconstructed signal ŝ. We observe a quite linear behavior. For 4 components the gain is of 17 db whereas for 16 components we gain 11.5 db. Finally, Fig. 5.6 (d) exposes the quadratic error of the peaks amplitude. There are two cases here. When all the components are found this error is linear (see Fig. 5.6 (d) the bottom curve when snr> 1) but it increases more quickly when some components are not found (see Fig. 5.6 (d) the bottom curve when snr< 1).

Conclusion.

In this paper, we have defined an iterative algorithm based on the descent gradient principle and adapted to the context of variational Bayesian methods. We also prove the convergence of this method in the probability density functions sets, ensuring that it converges to an optimal approximation of the posterior distribution. The main interest of this algorithm is that it converges faster than the classical Bayesian methods and allows a use on large dimensional datasets. We have furthermore give its implementation in the case of a white noise model when the prior information enhances some sparse behavior. A small tomographic application allows to compare our method with classical ones. We see that even in small cases, our algorithm can be faster than classical ones. A second simulation part, corresponding to a dictionary identification allows to understand the behavior of our method for large dimensional problems. Once again this method has good reconstruction properties in this case.

Appendix.

7.1. Proof of Lemma 2.3. In order to ensure that for small values of α, we have F (µ α ) ≥ F (µ k ), we show that the right part of Equation (2.17),

F (µ α ) -F (µ k ) ≥ ∂F µ k (µ α -µ k ) -L h α (µ k , .) -1 2 L 2 (µ k ) - R N ln(h α (µ k , x))dµ α (x),
is positive.

2.dF q k (q αq k ) -M 2 (q αq k ) 2 > 0, for all α < α 0 . This ensures that F (q α ) > F (q k ) for α ∈ (0, α 0 ).

First, one can notice that - Hence,

R N ln(h α (µ k , x))dµ α (x) = - R N (αdf (µ k , x) + ln K k (α))h α (µ k , x)dµ k (x) = -α R N df (µ k , x)h α (µ k , x)dµ k (x) + 1 α ln K k (α) .
∂F µ k (µ α -µ k ) - R N ln(h α (µ k , x))dµ α (x) ≥ R N (1 -α)df (µ k , x)(h α (µ k , x) -1)dµ k (x) = 1 -α α R N (αdf (µ k , x) + ln K k (α))(h α (µ k , x) -1)dµ k (x) - R N (1 -α) ln K k (α) α (h α (µ k , x) -1)dµ k (x) = 1 -α α R N
(αdf (µ k , x) + ln K k (α))(h α (µ k , x) -1)dµ k (x),

as (1-α) ln K k (α) α is constant and R N (h α (µ k , x) -1)dµ k (x) = 0. Finally, (αdf (µ k , x) + ln K k (α))( 1 α -1) -L(h α (µ k , x) -1) (h α (µ k , x) -1)dµ k (x).

∂F µ k (µ α -µ k ) -L h α -1 2 L 2 (µ k ) - R N
(7.3)

Let us consider each integrals appearing in Equation ( 7.3) separately. First, we can notice that if αdf (µ k , x) + ln K k (α) < 0, so is h α (µ k , .) -1. Furthermore, for every α > 0 and x ∈ R N , we have h α (µ k , .) -1 ≥ αdf (µ k , x) + ln K k (α). Hence if x ∈ R N is such that αdf (µ k , x) + ln K k (α) < 0 and α < 1 then

(αdf (µ k , x) + ln K k (α))( 1 α -1) -L(h α (µ k , x) -1) (h α (µ k , x)-1) ≥ (h α (µ k , x)-1) 2 1 α -1 -L ,
which is positive as soon as α ≤ 1 1+L .

Consider now that x ∈ R N is such that αdf (µ k , x) + ln K k (α) ≥ 0. The Mean Value Theorem applied to the exponential function ensures that one can find, for every x ∈ R N and α > 0, a θ(x, α) ∈ (0, αdf (µ k , x) + ln K k (α)) such that h α (µ k , x) = e αdf (µ k ,x)+ln K k (α) = 1+(αdf (µ k , x)+ln K k (α))e θ(x,α) ≥ 1+(αdf (µ k , x)+ln K k (α)).

This entails that

(αdf (µ k , x) + ln K k (α))( 1 α -1) -L(h α (µ k , x) -1) = (αdf (µ k , x)+ln K k (α))( 1 α -1-Le θ(x,α) ).

Furthermore, from Jensen's inequality 0 ≤ αdf (µ k , x) + ln K k (α) ≤ α df (µ k , x) -

R N df (µ k , x)dµ k (x) .
Thus, in this case,

df (µ k , x) - R N
df (µ k , y)dµ k (y) ≥ 0.

And

1 ≤ e θ(x,α) ≤ h α (µ k , x) ≤ e α(df (µ k ,x)-R N df (µ k ,x)dµ k (x)) ,

If we consider the function α → 1 α -1-Le α(df (µ k ,x)-R N df (µ k ,x)dµ k (x)) , it is decreasing and, as df is continuous, it is also positive for small values of α. Hence, one can find an α 1 > 0 such that for every α ≤ α 1 , 1 α -1 -Le α(df (µ k ,x)-R N df (µ k ,x)dµ k (x)) ≥ 0.

Finally one has

∀α ≤ α 0 = min(α 1 ,

1 1 + L ), ∂F µ k (µ α -µ k )-L h α -1 2 L 2 (µ k ) - R N
(αdf (µ k , x)+ln K k (α))dµ α (x) ≥ 0.

(7.4) 7.2. Context of simulation 5.2. In this part we present the dictionary decomposition of Section 5.2 together with the components of the data used to enhance the reconstruction performances of our unsupervised method.

To build our DFT-chirps dictionary we make the following assumptions : all chirps have the same duration T Chirp , the chirps rate ζ is digitalized on very few values (eight in the present case), and the Discrete Fourier Transform basis is included in the dictionary. In this context, the spurious signals can be represented on very few coefficients of the dictionary. However, we do not make any assumption on the variance of the noise or on the number of chirp functions in the mixture.

s(t) = N f req i=1 (u i + jv i )φ i (t) + Ntrans l=1 Nrate k=1 c l,k ψ k (t -l∆ t ), (7.5)
where

φ i (t) = exp [2jπf i t] ,
corresponds to a pure frequency f i with j 2 = -1 whereas

ψ k (t) = cos(2π(f 0 t + 1 2 t 2 ζ k ))Π 0,T Chirp (t)
corresponds to the chirps components. Here, f 0 is the starting frequency (at time t = 0), ζ k is the chirp rate, that is the increasing rate of the frequency, Π(t) is a gate function, T Chirp is the duration of the chirp, ∆ t is the shift between two chirps functions and t l = l∆ t is the first time where ψ k (t -l∆ t ) is not null. We merge all the dictionary's coefficients in a single vector x = (u 1 , . . . , u N f req , v 1 , . . . , v N f req , c 1,1 , . . . , c N trans,N rate ) t . Where N f req is the number of pure frequency functions contained on the dictionary whereas N trans is the number of chirp shifts and N rate is the number of chirp rates. The sampled version of Equation (7.5) is given by s = Hx, (

where s = (s(t 0 ), . . . , s(t N )) t . Observations are thus given by Equation (5.1).

The dictionary is composed of chirp functions with only 8 different chirp rates (N rate = 8). The frequency f 0 is 5 000 hz, the chirp rates ζ k are uniformly spaced between 6 000 and 20 000 hz, the shift parameter rate ∆ t is fixed up to a sampling period (T e = 1/F e ). Finally, the duration of the chirp (T Chirp ) is equal to the half time of the measurements (T mes ). Our dictionary is redundant as the number of coefficients is 4.5 times greater than the number of observations.

4 .

 4 Application to linear inverse problems. 4.1. Statement of the problem. The next part of this paper presents the application of Algorithm 2 to linear inverse ill-posed problems. The model of observations chosen in the following is given by y = Hx + b,

3 :function 4 : 5 :

 345 Estimate qk+1 (q k ) update a k+1 by Equation (4.7) update b k+1 by Equation (4.8)

  α ) using Equation (4.14) and Equation(4.13)

  .1 and Fig. 5.2(a)). Data have Coordinate (28,28) (25,28) (28,25) (40,28) (32,38) (48,48) (8

Fig. 5 . 1 .

 51 Fig. 5.1. Data collected : sinogram composed of 32 angles and 95 detector cells.

Fig. 5 . 2 .

 52 Fig. 5.2. Images are presented with the same inversed grayscale: (a) true image of 7 peaks, (b) FBP with ramp filter, (c) MCMC Gibbs approach, (d) classical variational Bayesian (VBBloc) with bloc optimization, (e) classical variational Bayesian (VBComp) with component optimization, (f ) SSR approach.

Fig. 5 . 3 .

 53 Fig. 5.3. Results with hyperparameters estimation: (a) True image, (b) reconstructed image with SSR algorithm (hyperparameters are fixed), (c) reconstructed image with USSR (image and hyperparameters are estimated jointly), (d) logarithm of σ 2 b : the dashed line correspond to the true value, (e) logarithm of σ 2 s

Fig. 5 . 4 .

 54 Fig. 5.4. Data in different fields: (a) time, (b) frequency, (c) time-frequency

Fig. 5 . 5 .

 55 Fig. 5.5. Limit case: (a) data and true signal, (b) reconstructed signal

Fig 5 . 6

 56 (a-d). Each point of the plot is computed by averaging the results of 30 reconstructions. In Fig 5.6 (a), resp. (b), we plot the true positive proportion, resp. false positive proportion, of the significant reconstructed coefficients * .

Fig. 5 . 6 .

 56 Fig. 5.6. Study of Robustness versus noise, each point of the curve being calculated by averaging 30 reconstructions with components randomly chosen: (a) True positive proportion of the significant reconstructed coefficients, (b) False positive proportion of the significant reconstructed coefficients, (c) snr of the estimated signal ŝ, (d) Quadratic error of the peaks amplitude

But

  Jensen's inequality ensures thatln K k (α) ≤ -α R N df (µ k , x)dµ k (x). α (µ k , x))dµ α (x) ≥ -α R N df (µ k , x)(h α (µ k , x) -1)dµ k (x).(7.2)

  ln(h α (µ k , x))dµ α (x) ≥ R N (αdf (µ k , x) + ln K k (α))( 1 α -1) -L(h α (µ k , x) -1) (h α (µ k , x) -1)dµ k (x) = {x:αdf (µ k ,x)+ln K k (α)>0} (αdf (µ k , x) + ln K k (α))( 1 α -1) -L(h α (µ k , x) -1) (h α (µ k , x) -1)dµ k (x) + {x:αdf (µ k ,x)+ln K k (α)≤0}

Table 5 . 2

 52 Comparison of the different approaches: computing time (second) and quality of estimation ( snr).

	Method	FBP VBBloc VBComp VBGrad (SSR) MCMC Gibbs
	CPU time (s) 0.05	586.20	1759.1	44.41	37079.50
	Nb of iterate	1	15	8(×4096)	500	1000
	snr	-2.04	5.87	5.76	6.00	-0.60
	by introducing a non-informative Jeffrey's prior, as described in Part 4.

Table 5 .

 5

			3
		Parameters of the different components
	type of function Amplitude frequency (hz) rate (hz) first time (s)
	cosine	1	5169.7	-	-
	cosine	0.8	4834	-	-
	Chirp	1.4	-	8000	0.2
	Chirp	1.4	-	10000	0.25
	Chirp	1.0	-	16000	0.22
	Chirp	1.0	-	20000	0.5
	Chirp	1.2	-	10000	0.4
	Chirp	1.0	-	18000	0.41
	Chirp	1.0	-	20000	0.6
	Chirp	1.4	-	8000	0.3

Table 5 .

 5

		4
		Parameters of the different components
	type of function Amplitude Chirp rate (hz) first time (s)
	Chirp	0.9526	16000	0.1530
	Chirp	1.1606	16000	0.1945
	Chirp	0.7369	18000	0.2000
	Chirp	1.1724	18000	0.1865

Table 5

 5

			.5
			Signal to Noise Ratio ( snr) in db
	data USSR approach Wiener filter wavelet shrinkage STFT shrinkage
	-5.0	15.05	1.1941	1.8677	9.2147

* the significant coefficients are obtained by applying a threshold equal to 0.2 on the coefficients vector. This threshold is equal to the third of the minimum value of the true non zero coefficients.