1. Introduction In many generic inverse problems in signal and image processing, the problem is to infer on an unknown signal f (t) or an unknown image f (r) with r = (x, y) through an observed signal g(t′ ) or an observed image g(r ′ ) related between them through an operator H such as convolution g = h ∗ f or any other linear or non linear transformation g = Hf . When this relation is linear and we have discretized the problem, we arrive to the relation: g = Hf + ǫ, where f = [f 1 , · · · , f n ]′ represents the unknowns, g = [g 1 , · · · , g m ]′ the observed data, ǫ = [ǫ1 , · · · , ǫm ]′ the errors of modelling and measurement and H the matrix of the system response. The Bayesian inference approach is based on the posterior law: p(f |g, θ 1 , θ 2 ) =

p(g|f , θ 1 ) p(f |θ 2 ) ∝ p(g|f , θ 1 ) p(f |θ 2 ) p(g|θ 1 , θ 2 )

(1)

where the sign ∝ stands for ”proportional to”, p(g|f , θ 1 ) is the likelihood, p(f |θ 2 ) the prior model, θ = (θ 1 , θ 2 ) are their corresponding parameters (often called the hyper-parameters of the problem) and p(g|θ 1 , θ 2 ) is called the evidence of the model. When the parameters θ have to be estimated too, a prior p(θ|θ 0 ) with fixed values for θ0 is assigned to them and the expression of the joint posterior p(g|f , θ 1 ) p(f |θ 2 ) p(θ|θ 0 ) p(f , θ|g, θ 0 ) = (2) p(g|θ 0 ) is used to infer them jointly. Variational Bayesian Approximation (BVA) methods try to approximate p(f , θ|g) by a e g) q2 (θ|f e , g) and then using them for estimation [1, 2, 3, separable one q(f , θ|g) = q1 (f |θ, 4, 5, 6, 7, 8, 9].

For hierarchical prior models with hidden variables z, the problem becomes more complex, because we have to give the expression of the joint posterior law p(f , z, θ|g) ∝ p(g|f , θ 1 ) p(f |z, θ 2 ) p(z|θ 3 ) p(θ|θ 0 )

(3)

and then approximate it by a separable one e g) q2 (z|f e , θ, e g) q3 (θ|z e , g) e , θ, e, f q(f , z, θ|g) = q1 (f |z

(4)

and then using them for estimation. In this paper, first the general VBA method is detailed for the inference on inverse problems with hierarchical prior models. Then, two particular classes of prior models (Student-t and mixture of Gaussians) are considered and the details of BVA algorithms for them are given. 2. Bayesian Variational Approximation with hierarchical prior models When a hierarchical prior model p(f |z, θ) is used and when the estimation of the hyperparameters θ has to be considered, the joint posterior law of all the unknowns becomes: p(f , z, θ|g) ∝ p(g|f , θ 1 ) p(f |z, θ 2 ) p(z|θ 3 ) p(θ)

(5)

The main idea behind the VBA is to approximate the joint posterior p(f , z, θ|g) by a separable one, for example q(f , z, θ|g) = q1 (f |g) q2 (z|g) q3 (θ|g) (6) and where the expressions of q(f , z, θ|g) is obtained by minimizing the Kullback-Leibler divergence Z q q KL(q : p) = q ln = ln (7) p p q It is then easy to show that KL(q : p) = ln p(g|M) − F(q) where p(g|M) is the likelihood of the model Z Z Z p(f , z, θ, g|M) df dz dθ (8) p(g|M) = with p(f , z, θ, g|M) = p(g|f , θ) p(f |z, θ) p(z|θ) p(θ) and F(q) is the free energy associated to q defined as p(f , z, θ, g|M) F(q) = ln (9) q(f , z, θ) q So, for a given model M, minimizing KL(q : p) is equivalent to maximizing F(q) and when optimized, F(q ∗ ) gives a lower bound for ln p(g|M). Without any other constraint than the normalization of q, an alternate optimization of F(q) with respect to q1 , q2 and q3 results in o n ) ∝ exp − hln p(f , z, θ, g)i q (f q(z )q(θ ) o , 1 n

q (z) ∝ exp − hln p(f , z, θ, g)i

2 q(f )q(θ ) o n q3 (θ) ∝ exp − hln p(f , z, θ, g)i q(f )q(z )

(10)

Note that these relations represent an implicit solution for q1 (f ), q2 (z) and q3 (θ) which need, at each iteration, the expression of the expectations in the right hand of exponentials. If p(g|f , z, θ 1 ) is a member of an exponential family and if all the priors p(f |z, θ 2 ), p(z|θ 3 ), p(θ 1 ), p(θ 2 ), and p(θ 3 ) are conjugate priors, then it is easy to see that these expressions leads

to standard distributions for which the required expectations are easily evaluated. In that case, we may note e g) q (z|f e , θ; e g) q (θ|f e, z e , θ; e ; g) q(f , z, θ|g) = q1 (f |z (11) 2 3

e and θ e are, respectively functions of (f e ,θ), e (z e and (f e ,z e, f e ,θ) e) where the tilded quantities z e e , θ) for and where the alternate optimization results to alternate updating of the parameters (z e e e e q1 , the parameters (f , θ) of q2 and the parameters (f , z ) of q3 . Finally, we may note that, to monitor the convergence of the algorithm, we may evaluate the free energy

F(q) = hln p(f , z, θ, g|M)iq + h− ln q(f , z, θ)iq = hln p(g|f , z, θ)iq + hln p(f |z, θ)iq + hln p(z|θ)iq + h− ln q(f )iq + h− ln q(z)iq + h− ln q(θ)iq (12) where all the expectations are with respect to q. Other decompositions are also possible: q(f , z, θ|g) =

Y j

or

Y

e g) e , θ; q1j (f j |fe (−j) , z

e g) e , θ; q(f , z, θ|g) = q1 (f |z

j

Y j

e g) e (−j), θ; q2j (zj |fe , z

e g) e (−j) , θ; q2j (zj |fe , z

Y l

Here, we consider this case and give some more details on it.

Y l

e e, θ q3l (θl |fe , z (−l) ; g)

e e, θ q3l (θ l |fe , z (−l) ; g)

3. Bayesian Variational Approximation with Student-t priors The Student-t model is: −(ν+1)/2 Y 1 Γ((ν + 1)/2) p(f |ν) = 1 + f 2j /ν St(f j |ν) with St(f j |ν) = √ πν Γ(ν/2) j Knowing that St(f j |ν) =

Z

0

(13)

(14)

(15)

∞

N (f j |0, 1/τ j ) G(τ j |ν/2, ν/2) dτ j

(16)

we can write this model via the positive hidden variables τ j : p(f |τ )

=

Q

j p(f j |τ j ) =

Q

(α−1)

p(τ j |α, β) = G(τ j |α, β) ∝ τ j

n

1 j N (f j |0, 1/τ j ) ∝ exp − 2

P

2 j τjfj

exp {−βτ j } with α = β = ν/2

o

(17)

Cauchy model is obtained when ν = 1: In this case, let consider the forward model g = Hf + ǫ and assign a Gaussian law to the noise ǫ which which results to p(g|f , vǫ ) = N (g|Hf , vǫ I). We also assign a prior also note τ =Q[τ1 , · · · , τN ], T = diag [τ ], p(τǫ |ατ 0 , βτ 0 ) = G(τǫ |ατ 0 , βτ 0 ) to τǫ = 1/vǫ . Let Q −1 zj = 1/τ j , Z = diag [z] = T and note p(f |τ ) = j p(f j |τ j ) = j N (f j |0, τ j ) = N (f |0, T ) Q and finally, assign p(τ |α0 , β0 ) = j G(τ j |α0 , β0 ). Then, we obtain the following expressions for the VBA:

p(g|f , τǫ ) = N (g|Hf , (1/τǫ )I) p(τǫ |ατ 0 , βQ τ 0 ) = G(τǫ |ατ 0 , βτ 0 ) p(f |τ ) = j N (f j |0, 1/τ j ) Q p(τ |α0 , β0 ) = j G(τ j |α0 , β0 ) e = N (f |µ, e e Σ) e Σ) q1 (f |µ, ′ e e = hλi ΣH g µ ′ e e −1 , Σ = (hλi H H + Z) e =T e −1 = diag [τe ] with Z

e j , βej ) q2j (τ j ) = G(τ j |α e j = α00 + 1/2 α e = β + < f 2 > /2 β 00 j j e τǫ , βeτǫ ), q3 (τǫ ) = G(τǫ |α α e τǫ = ατ 0 + (n + 1)/2 e β = βτ 0 + 1/2[kgk2 τǫ ′ ′ ′ ′

−2 hfi H g + H hf f i H]

e < f >= µ ′ e +µ eµ e′ < f f >= Σ 2 e e2j < f j >= [Σ]jj + µ e=α e τ /βeτ λ e j /βej τej = α

(18)

We can also express the free energy expression F(q) which can be used as a stopping criterion for the algorithm. Its expression is given in the appendix. The resulting algorithm can be summarized as follows e = N (f e , Σ) e e q1 (f |e τ , λ) λ −→ 6 e=λ eΣH e ′g f τe e e ′H + T e −1 )−1 −→ Σ = (λH 6

e αj , βej ) e q2j (τ j |f ) = G(τ j |e f −→ α ej = α00 + n+1

2 e βej = β00 + 1 f 2j Σ 2 −→ ej /βej τej = α

e f −→ e Σ −→ τej −→

e ) = G(τ |e q3 (τ |f ατ , βeτ ) n+1 α eτ = ατ 0 + 2 βeτ = βτ 0 + 12 [kgk2 − 2 < f >′ H ′ g + H ′ < f f ′ > H] e=α λ eτ /βeτ

e λ −→ τe −→

4. Bayesian Variational Approximation with Mixture of Gaussians priors The mixture models are also very commonly used as prior models. In particular the Mixture of two Gaussians (MoG2) model: p(f |λ, v1 , v0 ) =

Y j

(λN (f j |0, v1 ) + (1 − λ)N (f j |0, v0 ))

(19)

which can also be expressed through the binary valued hidden variables zj ∈ {0, 1}

p(f |z) =

Q

j

p(f j |zj ) =

P (zj = 1) = λ,

Q

j

N f j |0, vzj ∝ exp

P (zj = 0) = 1 − λ

− 21

P

f 2j j vzj

(20)

In general v1 >> v0 and λ measures the sparsity (0 < λ