domain: spi math

text:
Yuling Zheng
email: zheng@lss.supelec.fr

Aurélia Fraysse
email: fraysse@lss.supelec.fr

Thomas Rodet
email: trodet@satie.ens-cachan.fr

Efficient Variational Bayesian Approximation Method Based on Subspace optimization

Keywords: Variational Bayesian approximation, subspace optimization, large dimensional problem, unsupervised approach, total variation I

Variational Bayesian approximations have been widely used in fully Bayesian inference for approximating an intractable posterior distribution by a separable one. Nevertheless, the classical variational Bayesian approximation (VBA) method suffers from slow convergence to the approximate solution when tackling large-dimensional problems. To address this problem, we propose in this paper an improved VBA method. Actually, variational Bayesian issue can be seen as a convex functional optimization problem.

The proposed method is based on the adaptation of subspace optimization methods in Hilbert spaces to the function space involved, in order to solve this optimization problem in an iterative way. The aim is to determine an optimal direction at each iteration in order to get a more efficient method. We highlight the efficiency of our new VBA method and its application to image processing by considering an ill-posed linear inverse problem using a total variation prior. Comparisons with state of the art variational Bayesian methods through a numerical example show the notable improved computation time.

objects. The resolution of ill-posed inverse problems generally relies on regularizations which consist of introducing additional information, see [START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structure and problems[END_REF] for details. One most commonly used regularization is the Tikhonov one [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF]. Nevertheless, Tikhonov regularization leads to an estimator linear with respect to the data. Therefore, its ability to reconstruct non-linear components, such as location and magnitude of jumps or higher order discontinuities, see [START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF] and [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], is limited.

To overcome such limitations, nonlinear regularizations have been widely used. However, the drawback of those regularizations is the corresponding non quadratic or even non-convex optimization problems which are generally intricate. To tackle this issue, Geman et al [START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF], [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] proposed half-quadratic schemes in order to get nonlinear estimates more easily. By introducing auxiliary variables using duality tools, half-quadratic schemes transform the original complicated criterion into a half quadratic one where the original variables appear quadratically and the auxiliary variables are decoupled. This half-quadratic criterion can be efficiently optimized using classical optimization algorithms, which lead to the desired nonlinear estimates.

In a statistical framework, the half-quadratic schemes have been shown by Champagnat et al. [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF] as instances of the EM algorithm with latent variables which provide maximum a posteriori estimates.

Nevertheless, for either the Tikhonov regularization based methods or half-quadratic schemes, only point estimates could be given. Some useful information such as the variance of the estimator, which evaluates its precision, could not be directly obtained. However, such information is accessible if we obtain the posterior distribution of the unknown parameters, which is involved in the Bayesian framework.

Furthermore, another advantage of the Bayesian framework is that it provides a systematic way to determine hyperparameters, e.g. a Bayesian hierarchical approach can estimate the hyperparameters as well as unknown parameters by introducing hyperpriors, [START_REF] Molina | On the hierarchical Bayesian approach to image restoration: Application to astronomical images[END_REF] [START_REF] Higdon | Fully Bayesian estimation of Gibbs hyperparameters for emission computed tomography data[END_REF]. Such approaches are known as unsupervised approaches in the literature. In order to exploit these advantages, in the following we are mainly interested in the development of efficient unsupervised Bayesian reconstruction approaches.

Nevertheless, the difficulty met in general is that one could only acquire a posterior distribution whose partition function is unknown due to an intractable integral. In such a case, the main challenge is to retrieve the posterior distribution.

In this context, two main types of approaches are employed, stochastic approaches and analytic approximations. Stochastic approaches are based on Markov Chain Monte Carlo (MCMC) techniques [START_REF] Robert | Monte-Carlo Statistical Methods, ser. Springer Texts in Statistics[END_REF] which provide an asymptotically exact numerical approximation of the true posterior distribution.

The main drawbacks of such approaches are the high computational cost and the poor performance for large dimensional problems involving complicated covariance matrices. The use of such approaches for large dimensional problems is therefore limited.

Concerning analytic methods, MacKay in [START_REF] Mackay | Ensemble learning and evidence maximization[END_REF], see also [START_REF] Šmídl | The Variational Bayes Method in Signal Processing[END_REF] for a survey, proposed the variational Bayesian approximation (VBA) which aims to determine analytic approximations of the true posterior distribution. In this case, the objective is to find a simpler probability density function (pdf), generally separable, which is as close as possible to the true posterior distribution in the sense of minimizing the Kullback-Leibler divergence. This problem can be formulated as a convex infinite-dimensional optimization problem, whose resolution results in an optimal analytic approximation. However, this approximation does not have an explicit form except for extremely simple cases. In practice, it is generally approached by cyclic iterative methods which update at each iteration one component of the separable distribution while fixing the other ones. Such optimization procedure is known to be time consuming in general. The classical Bayesian methodology is thus not efficient enough when dealing with very large dimensional problems.

In order to obtain more efficient variational Bayesian approaches, a different method has been recently introduced in [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF]. It is based on the adaptation of the exponentiated gradient algorithm [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] into the space of pdf, which is no longer a Hilbert space. Instead of approximating an analytical solution of the involved functional optimization problem, this method seeks an approximate solution of this problem iteratively thanks to a gradient-type algorithm with explicit update equations. The optimization of the components of the separable approximation can thus be performed in parallel, which leads to a significant acceleration compared to the classical VBA method.

In order to further improve the method of [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF], a natural idea is to consider a new descent direction.

In this context, we propose to adapt the subspace optimization methods [START_REF] Miele | Study on a memory gradient method for the minimization of functions[END_REF]- [START_REF] Chouzenoux | A majorize-minimize subspace approach for l2-l0 image regularization[END_REF], into the space of pdf involved in variational Bayesian methodology. The advantage of subspace optimization is its generalized descent directions where Hilbert structure is not required. Moreover, the descent direction can be freely chosen in a subspace of dimension greater than one. This flexibility allows subspace optimization methods to be generally more efficient than conjugate gradient methods [START_REF] Chouzenoux | A Majorize-Minimize strategy for subspace optimization applied to image restoration[END_REF]. Based on the subspace optimization methods, an improved variational Bayesian optimization method is proposed in this paper. Moreover, we also consider the application of our improved variational Bayesian method to ill-posed linear inverse problems in image processing. In this context, the total variation (TV) regularization has been popular [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], [START_REF] Pollak | Nonlinear evolution equations as fast and exact solvers of estimation problems[END_REF] due to its ability to describe piecewise smooth images. Nevertheless, it is difficult to integrate the TV based prior into the development of unsupervised Bayesian approaches since its partition function depends on hyperparameters and it is not explicitly known. To tackle this problem, Bioucas-Dias et al. [START_REF] Bioucas-Dias | Adaptive total-variation image deconvolution: A majorization-minimization approach[END_REF], [START_REF] Oliveira | Adaptive total variation image deblurring: a majorizationminimization approach[END_REF] proposed a closed form approximation of this partition function. Thanks to this approximation, Babacan et al. [START_REF] Babacan | Variational Bayesian super resolution[END_REF] developed its TV based unsupervised Bayesian approach using the classical VBA method. In this work, we also take advantage of this approximate TV prior. With this prior, we develop our unsupervised approach using the proposed VBA method.

In order to evaluate the proposed approach, we provide also numerical comparisons with [START_REF] Babacan | Variational Bayesian super resolution[END_REF] which is based on the classical VBA method, and with another approach employing the gradient-type variational Bayesian algorithm proposed in [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF]. These comparisons are based on an implementation on a superresolution problem [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF] which aims to reconstruct a high resolution image from several low resolution ones representing the same scene. Moreover, in our reconstruction, we assume that motions between the low resolution images and a reference one could either be estimated in advance or be known through other sources of information. Such configuration appears for instance in astronomy [START_REF] Rodet | Data inversion for over-resolved spectral imaging in astronomy[END_REF] and medical imaging [START_REF] Greenspan | MRI inter-slice reconstruction using super-resolution[END_REF].

The rest of this paper is organized as follows: in Section II, we develop our proposed variational Bayesian optimization algorithm. Next, an application of the proposed algorithm to a linear inverse problem is shown in Section III whereas results of numerical experiments on super-resolution problems are given in Section IV. Finally, a conclusion is drawn in Section V.

II. EXPONENTIATED SUBSPACE-BASED VARIATIONAL BAYESIAN OPTIMIZATION ALGORITHM

A. Notations

In the following y ∈ R M and w ∈ R J denote respectively the data vector and the unknown parameter vector to be estimated whereas p(w), p(w|y) and q(w) represent the prior distribution, the true posterior law and the approximate posterior distribution, respectively.

B. Statement of the problem

The central idea of variational Bayesian methods is to approximate the true posterior distribution by a separable one

q(w) = P i=1 q i (w i ), (1)
where w = (w 1 , . . . , w P ). Here (w i ) i=1,...,P denote the P disjoint subsets of elements of w with P an integer between 1 and J.

The optimal approximation is determined by minimizing a measure of dissimilarity between the true posterior distribution and the approximate one. A natural choice for the dissimilarity measure is the Kullback-Leibler divergence (KL divergence), see [START_REF] Bernardo | Expected information as expected utility[END_REF]:

KL[q p(•|y)] = R J q(w) log q(w) p(w|y) dw. (2)
In fact, direct minimization of KL[q p(•|y)] is usually intractable since it depends on the true posterior p(•|y) whose normalization constant is difficult to calculate. However, as given in [START_REF] Choudrey | Variational methods for Bayesian independent component analysis[END_REF], the logarithm of the marginal probability of the data, called also evidence, can be written as

log p(y) = F(q) + KL[q p(•|y)], (3)
where F(q) is the so called negative free energy defined as

F(q) = R J q(w) log p(y, w) q(w) dw. (4)
As log p(y) is a constant with respect to q(w), minimizing the KL divergence is obviously equivalent to maximizing the negative free energy. We can see from (4) that for the computation of the negative free energy, the joint distribution p(y, w) is involved instead of the true posterior law. And this joint distribution can be readily obtained by the product of likelihood and prior distributions. We use hence the negative free energy as an alternative to the KL divergence.

Let us define a space Ω which is a space of separable pdfs, Ω = {q : pdf and q = P i=1 q i }. Our variational Bayesian problem is thus formulated as

q opt = arg max q∈Ω F(q) (5)
Classical variational Bayesian approximation [START_REF] Šmídl | The Variational Bayes Method in Signal Processing[END_REF] is based on an analytical solution of (5) which is

given by q = P i=1 q i with q i (w i ) = K i exp log p(y, w) j =i qj(wj) , ∀i = 1, . . . , P.

Here • q = E q [•] and K i denotes the normalization constant. We can see from ( 6) that q i depends on the other marginal distributions q j for j = i. As a result, we cannot obtain an explicit form for q unless in extremely simple cases. Therefore, iterative methods such as the Gauss-Seidel one have to be used to iteratively approach this solution. As a result, classical variational Bayesian method is not efficient enough to treat high dimensional problems.

In this paper, in order to obtain efficient variational Bayesian approaches, instead of firstly giving an analytical solution then iteratively approaching it, we directly propose iterative methods to solve [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF] which is a functional optimization problem over functions (q i ) i=1,...,P , in the space Ω. As stated in [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF],

there exists a problem equivalent to (5) in a separable probability measure space A = P i=1 A i , the Cartesian product of the A i , which is defined as follows:

A i = {µ i : probability measure and µ i (dw i ) = q i (w i )dw i with q i a pdf}.

Therefore, the optimization of ( 5) is equivalent to the resolution of the following optimization problem

µ opt = arg max µ∈A F (µ), (7)
where the functional F satisfies that ∀µ ∈ A of density q, F (µ) = F(q).

In [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF], the gradient descent method in Hilbert spaces has been transposed into the space of pdfs and, as a result, an exponentiated gradient based variational Bayesian approximation (EGrad-VBA) method whose convergence is proven was proposed to solve the involved functional optimization problem. For the aim of developing more efficient methods, we transpose here in the same context the subspace optimization method which has been shown to outperform standard optimization methods, such as gradient or conjugate gradient methods, in terms of rate of convergence in finite dimensional Hilbert spaces [START_REF] Chouzenoux | A Majorize-Minimize strategy for subspace optimization applied to image restoration[END_REF].

C. Subspace optimization method in Hilbert spaces

We give in this section a brief introduction of the subspace optimization method in Hilbert spaces. The subspace optimization method has been proposed by Miele et al. [START_REF] Miele | Study on a memory gradient method for the minimization of functions[END_REF] with a subspace spanned by the opposite gradient and the previous direction. This method is known as Memory Gradient (MG) and can be regarded as a generalization of the conjugate gradient method. More recently, a lot of other subspace optimization methods based on different subspaces, see [START_REF] Narkiss | Sequential Subspace Optimization Method for Large-Scale Unconstrained Problems[END_REF] and [START_REF] Shi | A new super-memory gradient method with curve search rule[END_REF] for instance, have been proposed.

Generally, subspace optimization methods use the following iteration formula:

x k+1 = x k + d k = x k + D k s k , (8)
where x k and x k+1 respectively stands for the estimates at kth and (k+1)th iterations, An overview of existing subspace optimization methods [START_REF] Chouzenoux | A Majorize-Minimize strategy for subspace optimization applied to image restoration[END_REF] shows that D k usually includes a descent direction (e.g. gradient, Newton, truncated Newton direction) and a short history of previous directions.

D k = [d k 1 , . . . , d k I ] gathers I directions
In this work, we consider only the super memory gradient subspace where D k is given as follows:

D k = [-g k , d k-1 , . . . , d k-I+1 ]. (9)
Here -g k is the opposite gradient and (d k-j ) j=1,...,I-1 are the directions of previous iterations. Chouzenoux et al. [START_REF] Chouzenoux | A Majorize-Minimize strategy for subspace optimization applied to image restoration[END_REF] have addressed a discussion about the dimension of the subspace through simulation results on several image restoration problems. It is shown that in a Hilbert space, for a super memory gradient subspace (9), taking I = 2 i.e. a subspace constructed by the opposite gradient and the previous direction, results in the best performance in terms of computation time. In this case, the super memory gradient subspace is degraded to the memory gradient one.

D. Proposed subspace-based variational Bayesian approximation method

In this section, we define our iterative method based on the transposition of the subspace optimization method for the resolution of [START_REF] Higdon | Fully Bayesian estimation of Gibbs hyperparameters for emission computed tomography data[END_REF]. We use here k ∈ N, set initially to zero, as the iteration count and assume that µ k is a Radon probability measure [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF] with a density q k , i.e. µ k (dw) = q k (w)dw. As we stand in the space of probability density measures, the next iteration should give also a probability density measure absolutely continuous with respect to µ k . The Radon-Nikodym theorem [START_REF] Rudin | Real and Complex Analysis[END_REF] ensures that this measure should be written as

µ k+1 (dw) = h k (w)µ k (dw), (10)
where

h k ∈ L 1 (µ k ) is a positive function 1 .
Since µ k is a Radon probability measure with a density q k , we can equivalently write

q k+1 (w) = h k (w)q k (w), (11)
as updating equation for the approximate density. Furthermore, as we deal with entropy-type functionals, a natural choice for h k would be an exponential form (see [START_REF] Tseng | On the convergence of the exponential multiplier method for convex programming[END_REF]). Moreover, this choice ensures the positivity of pdfs along iterations.

Considering the exponential form of h k and the subspace optimization principle, we propose

h k (w) = K k (s k ) exp D k (w)s k , (12)
where K k (s k ) represents the normalization constant expressed as

K k (s k ) = R J exp D k (w)s k q k (w)dw -1 , (13)
and

D k (w) = [d k 1 (w), . . . , d k I (w)]
is the set of I directions spanning the subspace. We should state that as we deal with a functional optimization problem, the directions (d k l (w)) l=1,...,I , are no longer vectors but functions. Finally, s k = [s k 1 , . . . , s k I ] T denotes the multi-dimensional step-size.

Due to the exponential form, [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] can also be written as:

h k (w) = K k (s k ) φ k 1 (w) s k 1 . . . φ k I (w) s k I ( 14
)
where φ k l (w) = exp[d k l (w)], for l = 1, . . . , I. 1) Set of directions spanning the subspace: As discussed in Section II-C, for the super memory gradient subspaces defined in Hilbert spaces (see [START_REF] Mackay | Ensemble learning and evidence maximization[END_REF]), the subspace of dimension two, which is known as memory gradient subspace, leads to the most efficient approaches. Moreover, a subspace of dimension two would result in less computation complexity than higher order subspaces. As a result, we consider in this work the transposition of the memory gradient subspace into the space of pdfs. In this case, I = 2.

The transposition leads to a D k (w) consisting of one term related to the Gateaux differential of F(q) and the other term corresponding to the previous direction. The Gateaux differential of the negative free energy F(q) is a linear functional defined on Ω by

∀(q, q) ∈ Ω 2 dF q (q) = R J df (q, w)q(w)dw. ( 15
)
In the case of separable q, we have df (q, w) = Σ i d i f (q, w i ), in which ∀i = 1, . . . , P , d i f (q, w i ) is expressed as:

d i f (q, w i ) = log p(y, w) j =i qj(wj) -log q i (w i ) -1. (16)
Mathematically, the structure of our Memory Gradient (MG) subspace is given by

D k M G (w) = [df (q k , w), d k-1 (w)], (17)
where df (q k , •) is given by ( 16) and d k-1 stands for the direction of the previous iteration, which is given by

d k-1 (w) = log i q k i (w i ) i q k-1 i (w i ) . (18)
Our proposed variational Bayesian approximation method based on this subspace is called exponentiated memory gradient subspace-based variational Bayesian approximation in the rest of this paper and is resumed as follows:

Concerning the step size, let us first define

g k : R 2 → R g k (s) = F K k (s)q k (w) exp D k M G (w)s , (19)
then the optimal step-size is given by

(ŝ opt ) k = arg max s∈R 2 g k (s). ( 20
)
Algorithm Exponentiated Memory Gradient subspace-based Variational Bayesian Approximation (EMG-VBA)

1) Initialize(q 0 ∈ Ω)

2) repeat a. determine subspace D k M G (w) using ( 16), ( 17) and ( 18) b. determine step-sizes s k c. compute

q k+1 (w) = K k s k q k (w) exp D k M G (w)s k until convergence
For the proposed EMG-VBA algorithm, we would prove the following proposition.

Proposition 1. Let q 0 ∈ Ω and ∀k ≥ 0, let the sequence {q k } k≥0 be defined by

q k+1 (w) =K k s k q k (w) exp D k M G (w)s k =K k s k q k (w) × exp s k 1 df (q k , w) + s k 2 d k-1 (w) . ( 21
)
If s k is the optimal step-size defined in [START_REF] Bioucas-Dias | Adaptive total-variation image deconvolution: A majorization-minimization approach[END_REF] then

{F(q k )} k≥0 converges to a maximum of F(q).
In a previous work [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF], it has been proven that for q k ∈ Ω, if q k+1 is constructed by the EGrad-VBA with the following updating formula:

q k+1 grad (w) = K k s k 1 q k (w) exp s k 1 df (q k , w) , (22)
with an optimal step-size defined as:

s k 1 = arg max s1∈R F K k (s 1 ) q k (w) exp s 1 df (q k , w) , (23)
then the negative free energy F(q) increases. More precisely,

F(q k+1 grad ) ≥ F(q k )
for q k ∈ Ω and q k+1 grad constructed by [START_REF] Babacan | Variational Bayesian super resolution[END_REF] (24)

Comparing ( 21) and ( 22), we notice that the EGrad-VBA algorithm can be identified as a special case of our proposed EMG-VBA algorithm with the second step-size s 2 set to zero. Due to the use of optimal May 12, 2014 DRAFT step-sizes, we obtain

F(q k+1 grad ) = max s1∈R F K k (s 1 ) q k (w) exp s 1 df (q k , w) = max s1∈R,s2=0 g k (s) (25)
However, in our proposed EMG-VBA, we consider

F(q k+1 ) = max s∈R 2 g k (s) ≥ max s1∈R,s2=0 g k (s) = F(q k+1 grad ) (26)
Combining ( 24) and ( 26), we obtain

F(q k+1 ) ≥ F(q k ), (27)
which shows that for a sequence {q k } k≥0 constructed by the proposed EMG-VBA, the sequence of negative free energy {F(q k )} k≥0 is a increasing one. Furthermore, F is a concave functional. As a result, the proposed EMG-VBA generates a sequence {F(q k )} k≥0 which converges to a maximum of

F(q).

E. The approximate step-size

Generally, it is too expensive to get the optimal step-size. Therefore, in practice most iterative approaches take sub-optimal ones considering the trade-off between computational cost and difference from the optimal step-size. In scalar cases, typical line search approaches determine trade-off steps by trying out a sequence of values until the fulfillment of certain sufficient conditions, such as Wolfe,

Goldstein [START_REF] Nocedal | Numerical Optimization, ser. Series in Operations Research[END_REF]. An extension of these conditions to multi-dimensional cases can be easily obtained.

However, such approximate methods in scalar cases have already been known to be time-consuming when the computation of the objective criterion is expensive, which is the case here. Moreover, it is difficult to adjust multi-dimensional step-sizes to satisfy the sufficient conditions. As a result, the extension to multi-dimensional cases can greatly increase the computational cost. Furthermore, the rate of convergence of such line search methods depends closely on parameters controlling the boundaries of chosen conditions and on the starting point of the step-size. A bad choice of these parameters will cause a slow convergence. Therefore, in this paper, we do not consider such approximate methods and take sub-optimal steps calculated as in [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF]. Firstly, we take the second order Taylor expansion of g k (s)

at origin as local approximation,

gk (s) = g k (0) + ∂g k ∂s s=0 T s + 1 2 s T ∂ 2 g k ∂s∂s T s=0 s, (28)
where ∂g k ∂s s=0 denotes the gradient vector whereas ∂ 2 g k ∂s∂s T s=0 is the Hessian matrix of the function g k at s = 0. When s 1 and s 2 are small, gk (s) is a close approximation of g k (s). Secondly, we compute our sub-optimal steps by maximizing gk (s) which is quadratic. Assuming that the Hessian matrix ∂ 2 g k ∂s∂s T s=0 is invertible, our sub-optimal steps are given by

(ŝ subopt ) k = - ∂ 2 g k ∂s∂s T s=0 -1 ∂g k ∂s s=0 . (29)

III. APPLICATION TO A LINEAR INVERSE PROBLEM WITH A TOTAL VARIATION PRIOR

We show in this section an application of the proposed EMG-VBA to ill-posed inverse problems in image processing. Variational Bayesian approaches are widely used to tackle inverse problems where complicated posterior distributions are involved. In the following, we firstly present such a problem which adopts a Total Variation (TV) prior. Then we develop an unsupervised Bayesian approach with the proposed EMG-VBA.

A. Direct model

We consider here a classical linear direct model:

y = Ax + n, (30)
where y ∈ R M and x ∈ R N denote respectively data and unknown parameters to be estimated arranged in column lexicographic ordering. The linear operator A ∈ R M ×N is assumed to be known and n is an additive white noise, assumed to be i.i.d. Gaussian, n ∼ N (0, γ -1 n I), with γ n as a precision parameter, i.e. the inverse of the noise variance. The direct model [START_REF] Tseng | On the convergence of the exponential multiplier method for convex programming[END_REF] and the hypothesis of i.i.d. Gaussian noise allow an easy derivation of the likelihood function:

p(y|x, γ n ) ∝ γ M/2 n exp - γ n y -Ax 2 2 . (31)

B. Image model

In this work, we consider an image model, more precisely a prior distribution of the unknown x, satisfying two main properties. Firstly, it is able to describe the piecewise smoothness property of images.

In the literature, total variation has been largely used in various imaging problems including denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], blind deconvolution [START_REF] Chan | Total variation blind deconvolution[END_REF], inpainting [START_REF] Shen | Mathematical models for local nontexture inpaintings[END_REF] and super-resolution [START_REF] Babacan | Variational Bayesian super resolution[END_REF]. Secondly, we should have some knowledge of its partition function in order to develop unsupervised approaches which sidestep the difficulty of tuning hyperparameters. Both demands described above lead us to the work of Babacan et al. [START_REF] Babacan | Variational Bayesian super resolution[END_REF], where an unsupervised Bayesian approach using the total variation (TV) prior was developed.

The TV prior is given by

p(x|γ p ) = 1 Z T V (γ p ) exp [-γ p T V (x)] , (32)
where Z T V (γ p ) is the partition function and

T V (x) = N i=1 (D h x) 2 i + (D v x) 2 i . (33)
Here D h and D v represent first-order finite difference matrices in horizontal and vertical directions, respectively.

The major difficulty is that there is no closed form expression for the partition function Z T V (γ p ). To overcome this difficulty, Bioucas-Dias et al. [START_REF] Bioucas-Dias | Adaptive total-variation image deconvolution: A majorization-minimization approach[END_REF], [START_REF] Oliveira | Adaptive total variation image deblurring: a majorizationminimization approach[END_REF] proposed an approximate partition function

Z T V (γ p ) ≃ Cγ -θN p , ( 34
)
where C is a constant and θ is a parameter which has to be adjusted in practice to get better results.

This analytic approximation results in

p(x|γ p ) ≃ p(x|γ p ) = Cγ θN p exp [-γ p T V (x)] . (35)
Babacan et al. [START_REF] Babacan | Variational Bayesian super resolution[END_REF] adopted this approximate TV prior with θ = 1/2. In such a case, Z T V (γ p ) is approximated by Cγ -N/2 p which corresponds to the partition function of a multivariate normal distribution of a N-dimensional random vector. However, a TV prior is not similar to a Gaussian one. Therefore, this approximation is not close enough. As a result, in this paper, we keep the parameter θ and adjust it in practice.

C. Hyperpriors

The hyperparameters γ n and γ p play an important role in the performance of algorithms. In practice, choosing correct hyperparameters is far from a trivial task. Therefore, we prefer to automatically determine their values. This is done by introducing hyperpriors for the hyperparameters. In order to obtain numerically implementable approaches, conjugate hyperpriors are employed. For γ n and γ p , we use Gamma distributions,

p(γ n ) = G(γ n |ã n , bn ) (36) p(γ p ) = G(γ p |ã p , bp ) (37)
where for a > 0 and b > 0

G(z|a, b) = b a Γ(a) z a-1 exp (-bz) (38)
As we do not have any prior information about γ n and γ p , in practice, we consider ãn = 0, bn = 0 and ãp = 0, bp = 0, which lead to non-informative Jeffreys' priors.

Consequently, we obtain a joint distribution as follows

p(y,x, γ n , γ p ) ∝ γ M/2 n exp - γ n y -Ax 2 2 × γ θN p exp -γ p N i=1 (D h x) 2 i + (D v x) 2 i γ -1 n γ -1 p (39)
where ∝ means "is approximately proportional to". The posterior distribution p(x, γ n , γ p |y) is not known explicitly since the normalization constant is not calculable. In order to proceed the statistic inference of the unknown variables, we resort to the variational Bayesian approximation methods which aims at getting the best separable analytical approximation. In the context of variational Bayesian methods, in order to get numerically implementable approaches, conjugate priors are needed to ensure that each posterior distribution belongs to a given family. Consequently, the optimization of the approximate distribution can be reduced to a numerical approximation of its parameters. Nevertheless, the TV prior introduced above is not conjugate with the likelihood (see [START_REF] Nocedal | Numerical Optimization, ser. Series in Operations Research[END_REF]) which is a Gaussian distribution. To tackle this problem, Minorization-Maximization (MM) techniques [START_REF] Hunter | A tutorial on MM algorithms[END_REF] are employed here to get a conjugate variant, as done by Babacan et al. in [START_REF] Babacan | Variational Bayesian super resolution[END_REF].

D. Application of variational Bayesian approximation methods

Let Θ denote the set of all unknown parameters: Θ = {x, γ n , γ p }, the objective of variational Bayesian approximations is to give a tractable approximation q Θ of the true posterior distribution p(•|y). As the TV prior is not conjugate with the likelihood, it is difficult to carry out the maximization of the negative free energy with respect to q Θ . This difficulty has been solved by adopting Minorization-Maximization (MM) techniques [START_REF] Hunter | A tutorial on MM algorithms[END_REF], in which maximizing the negative free energy is substituted by maximizing a tractable lower bound. To get such a lower bound of the negative free energy, a lower bound of the approximate TV prior is firstly constructed by introducing positive auxiliary variables λ = [λ 1 , . . . , λ N ], see [START_REF] Babacan | Variational Bayesian super resolution[END_REF] for details,

p(x|γ p ) ≥M (x, γ p |λ) = cγ θN p × exp -γ p N i=1 (D h x) 2 i + (D v x) 2 i + λ i 2 √ λ i . ( 40
)
From (40), we can see that the lower bound of the TV prior, M (x, γ p |λ), is proportional to a Gaussian distribution and is therefore conjugate to the likelihood. Combining ( 4) and ( 40), a lower bound of the negative free energy can be derived as

F(q Θ ) ≥ F L (q Θ , λ) = q Θ (x, γ n , γ p ) log L(x, γ p , γ n , y|λ) q Θ (x, γ n , γ p ) dxdγ n dγ p , (41)
where

L(x, γ p , γ n , y|λ) = p(y|x, γ n )M (x, γ p |λ)p(γ n )p(γ p ) (42)
is a lower bound of the joint distribution.

Hence the resolution of the problem ( 5) is performed by alternating the two following steps: maximizing the lower bound F L with respect to the pdf q Θ and updating the auxiliary variable λ in order to maximize

F L . Assuming that q Θ (Θ) = q x (x)q γn (γ n )q γp (γ p ) = i q i (x i )q γn (γ n )q γp (γ p ), (43)
we carry out an alternate optimization of F L with respect to q x , q γn , q γp and λ. Altogether, we perform the following alternate iterative scheme 2 :

q k+1 x = arg max qx F L q x q k γn q k γp , λ k , ( 44
)
λ k+1 = arg max λ∈R N F L q k+1 x q k γn q k γp , λ (45)
q k+1 γn = arg max

qγ n F L q k+1 x q γn q k γp , λ k+1 , (46)
q k+1 γp = arg max qγ p F L q k+1 x q k+1 γn q γp , λ k+1 (47)
The functional optimizations with respect to q x , q γn and q γp (given by ( 44), ( 46) and (47), respectively) are solved by variational Bayesian approximation methods. Since the conditional posterior p(γ n , γ p |x, y) is separable, it could be approximated efficiently thanks to the classical VBA. In fact, the proposed EMG-VBA is only adopted to approximate the posterior distribution of x where it improves the rate of convergence. As regards the optimization of the auxiliary variable λ (given by ( 45)), it involves a classical optimization in a Hilbert space. 2 The auxiliary variable λ is chosen to be updated before qγ n and qγ p in order to get a simpler iteration formula for qγ p , see [START_REF] Babacan | Variational Bayesian super resolution[END_REF].

May 12, 2014 DRAFT Due to the use of MM techniques, we manage to get a prior for x conjugate with the likelihood.

Moreover, conjugate Gamma priors are chosen for hyperparameters γ n and γ p . Therefore, the optimal approximations (q i ) i=1,...,N belong to a Gaussian family whereas the optimal approximate posterior distributions of hyperparameters q γn and q γp belong to a Gamma one.

q k x (x) = i N (x i |(m k ) i , (σ 2 k ) i ), (48)
q k γn (γ n ) = G(γ n |a k γn , b k γn ), (49)
q k γp (γ p ) = G(γ p |a k γp , b k γp ), (50)
Therefore, the optimization of approximate distributions can be performed by iteratively updating their parameters.

1) Optimization of q x using the proposed EMG-VBA: According to ( 11), ( 12) and ( 16) -( 18), we get a distribution of x depending on the step-size s:

q s x (x) = K k (s)q k x (x) exp(s 1 df (q k x , x) + s 2 d k-1 (x)) = K k (s)q k x (x) i q r i (x i ) q k i (x i ) s1 q k i (x i ) q k-1 i (x i ) s2 . ( 51
)
As the lower bound of the negative free energy is involved, the auxiliary function q r i is of the form

q r i (x i ) ∝ exp log L(x, γ p , γ n , y|λ k ) ( j =i q k j )q k γn q k γp ∝ exp - γ n 2 y -Ax 2 + γ p N i=1 (D h x) 2 i + (D v x) 2 i + λ k i 2 λ k i × j =i q k j (x j )dx j q k γn (γ n )q k γp (γ p )dγ n dγ p ∝ exp - γ n k 2 x 2 i diag A T A i -2x i A T y i + 2x i A T Am k i -2x i diag A T A i (m k ) i - γ p k 2 x 2 i diag D T h Λ k D h + D T v Λ k D v i + 2x i D T h Λ k D h m k + D T v Λ k D v m k i -2x i diag D T h Λ k D h + D T v Λ k D v i (m k ) i (52)
where

z k = E q k z (z), Λ k = Diag 1 √ λ k i is a diagonal matrix with 1 √ λ k i i=1,...,N
as its elements.

Moreover, diag (M) is a vector whose entries are the diagonal elements of matrix M.

The computation of q r i , for i = 1, . . . , N , shows that each of them can be identified as a Gaussian distribution with mean and variance expressed explicitly by the two following expressions:

(σ 2 r ) i = γ n k diag A T A i + γ p k diag D T h Λ k D h + D T v Λ k D v i -1 , (53)
(m r ) i =(σ 2 r ) i γ n k A T y -A T Am k + diag(A T A) • m k i -γ p k D T h Λ k D h m k + D T v Λ k D v m k i + γ p k diag D T h Λ k D h + D T v Λ k D v i (m k ) i , (54)
where • denotes the Hadamard product between two vectors.

Based on the above results for q r i , using (51), we can derive the expression of q s x (x) = i q s i (x i ) where each component q s i (x i ) is a Gaussian distribution with mean (m s ) i and variance (σ 2 s ) i satisfying:

σ 2 s = 1 σ 2 k + s 1 1 σ 2 r - 1 σ 2 k + s 2 1 σ 2 k - 1 σ 2 k-1 -1 , (55)
m s = σ 2 s m k σ 2 k + s 1 m r σ 2 r - m k σ 2 k + s 2 m k σ 2 k - m k-1 σ 2 k-1 . (56)
In above equations, we omit all the indication of vector component (•) i for the sake of clarity. From (55), we can see that σ 2 s is equal to the inverse of a linear combination of three terms where the first term is the present inverse variance, the second term comes from the gradient and the third term is caused by the memory of the previous direction. From (56), we can see that m s has the same structure. As stated earlier, the EGrad-VBA can be identified as the proposed EMG-VBA with s 2 set to zero which leads to the elimination of the third term in (55) and ( 56). Because of the extra term (the third term) compared to EGrad-VBA, EMG-VBA can obtain a closer approximation than EGrad-VBA in one iteration.

The previous distribution is still a function of the step size. A sub-optimal step-size defined by [START_REF] Rudin | Real and Complex Analysis[END_REF] in Section II-E is then adopted. As a result, (σ 2 ) k+1 = σ 2 ŝsubopt and m k+1 = m ŝsubopt . 2) optimization of λ in Hilbert spaces: The elements of auxiliary vector λ are calculated by maximizing the upper bound F L with respect to (λ i ) i=1,...,N . Since F L is concave and differentiable with May 12, 2014 DRAFT respect to (λ i ) i=1,...,N , its maximum is achieved at the critical point which is given by

λ k+1 i =E q k+1 x (D h x) 2 i + (D v x) 2 i =E q k+1 x x T (D h ) T i (D h ) i x + x T (D v ) T i (D v ) i x =(D h m k+1 ) 2 i + (D v m k+1 ) 2 i + trace (D h ) T i (D h ) i Σ k+1 ) + trace (D v ) T i (D v ) i Σ k+1 ) , (57)
where (D h ) i and (D v ) i represent the ith row of D h and D v , respectively. And

Σ k+1 = Diag (σ 2 ) k+1
is the covariance matrix which is diagonal under the separability hypothesis.

3) optimization of q γn and q γp using classical VBA: These two distributions are computed using [START_REF] Molina | On the hierarchical Bayesian approach to image restoration: Application to astronomical images[END_REF].

More details of the calculus can be found in [START_REF] Babacan | Variational Bayesian super resolution[END_REF]. The means of Gamma distributions are used as the estimates of hyperparameters which are

γ n k+1 = M E q k+1 x [ y -Ax 2 ] , (58)
γ p k+1 = θN N i=1 λ k+1 i . (59)
Altogether, the proposed algorithm for inverse problem using TV prior is summed up in Algorithm 1.

Algorithm 1 Proposed unsupervised variational Bayesian approach 1) Initialize parameters of (q 0 i ) i=1,...,N , q 0 γn , q 0 γp and λ

2) Update means and variances of q k+1 i for i = 1, . . . , N a. Compute parameters of intermediary functions q r i using ( 53 ) using ( 29) c. Update means and variances of q k+1 i using (55), (56)

3) Update auxiliary vector λ k+1 using (57) 4) Determine the parameters of q k+1 γn and compute its mean by (58) 5) Determine the parameters of q k+1 γp and compute its mean by (59) 6) Go back to 2) until convergence

IV. EXPERIMENTAL EVALUATION

The performance of the proposed approach (Algorithm 1) is evaluated through an application on a super-resolution (SR) problem. A SR problem is covered by the linear direct model (30) with a system matrix A gathering the warping, blurring and down-sampling operators. In fact, the main concern of this section is the evaluation of the time efficiency of the new variational Bayesian algorithm EMG-VBA by comparisons with the existing variational Bayesian approximation techniques, classical VBA and EGrad-VBA [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF] for the estimation of images. Therefore, we treat here a non-classical super-resolution problem where, for the sake of simplicity, the system matrix A is assumed to be known, i.e. no motion parameters are estimated. Meanwhile, this assumption would reduce the implementation limitation of a state of the art approach based on classical VBA [START_REF] Babacan | Variational Bayesian super resolution[END_REF] to large dimensional problems. In the following, at first, we present briefly this state of the art approach based on the classical VBA and another one using the EGrad-VBA, then we show the comparisons of the proposed approach with these two approaches.

A. State of the art approaches

1) A SR approach based on classical VBA [START_REF] Babacan | Variational Bayesian super resolution[END_REF]: The linear inverse problem with TV image prior treated in Section III has also been treated recently by Babacan et al. in [START_REF] Babacan | Variational Bayesian super resolution[END_REF] in the context of classical super-resolution. In this paper, we suppose that the motion parameters are known in order to make the SR approach of Babacan et al. applicable for large-dimensional problems.

Two major differences exist between the work in [START_REF] Babacan | Variational Bayesian super resolution[END_REF] and our work presented in Section III-D. The first one is that Babacan et al. used the classical VBA for the optimization of q x whereas we adopt the proposed EMG-VBA. The second difference is that we assume that q x is fully separable whereas Babacan et al. supposed that it is not. Due to the non-separability assumption, classical VBA based on (6) yields a multivariate Gaussian distribution for q k x . Updating the distribution is also equivalent to updating its parameters, i.e. mean m k and covariance matrix Σ k . These two parameters are given by

m k+1 =Σ k+1 γ n k A T y , (60)
(Σ k+1 ) -1 = γ n k A T A + γ p k D T h Λ k D h + D T v Λ k D v , (61)
We can see from (60) that m k+1 depends on the covariance matrix Σ k+1 but the computation of Σ k+1 needs the inversion of the matrix given by (61). To bypass the matrix inversion, Babacan et al. have adopted the conjugate gradient method to iteratively approximate m k+1 , which can be inefficient in large dimensional cases. For the optimization of the auxiliary variable λ and the hyperparameters γ n , γ p , the same updating equations as (57), ( 58) and (59) are obtained. In (57) where the covariance matrix Σ k+1 is needed, it is approximated by a diagonal matrix whose diagonal entries are equal to the inverse of the diagonal elements of (Σ k+1 ) -1 . Generally, this approximation is not a precise one.

2) EGrad-VBA [START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF] based approach: We present in this section another approach for comparison which treats the same inverse problem. The only major difference between this approach and the proposed one is that the EGrad-VBA using approximate optimal step-sizes defined in a same way as ( 29) is adopted for the optimization of q x . In fact, EMG-VBA differs from the EGrad-VBA only in the selection of the direction: EMG-VBA takes the memory of previous direction into consideration whereas EGrad-VBA does not. Here, there is no need to show the details of the application of EGrad-VBA into our inverse problem since it is a degenerated version of the EMG-VBA with a subspace D k (w) consisting of only df (q k , •). Consequently, the updating equations of the mean and variance of q x could be easily obtained by considering s 2 = 0 in (55) and (56). The updating equations for the auxiliary variable λ and the hyperparameters γ n , γ p are still the same as (57), ( 58) and (59).

B. Simulation results

In the following, the approach of Babacan et al. [START_REF] Babacan | Variational Bayesian super resolution[END_REF] dB. However, for the sake of simplicity, we treat a problem where LR images are without rotation and motions are supposed to be known. During the reconstruction, we take twelve LR images as data and assume that the convolution kernel, decimation factor and shifting of LR images are all known. Since the decimation factor of 4 is used, the size of a LR image is 1 16 of that of the HR image. As a result, the size of twelve LR images is smaller than that of the objective HR image.

In fact, all the considered approaches are based on a same Bayesian model and tackle the same optimization problem. Therefore, in general, these approaches lead to similar reconstruction results. The main concern of this comparison is their rate of convergence. To have a fair comparison, we take the same initializations for all the approaches: m 0 = A T y as the mean and 100 as the variance of HR May 12, 2014 DRAFT image pixels, the initializations of the auxiliary variables and the hyperparameters are computed from m 0 using (57), ( 58) and (59). Moreover, the parameter θ involved in the partition function of the image prior needs to be adjusted. A set of experiments carried out with different images show that the best results are achieved with θ = 1.1. As a result, we set θ = 1.1 for all our experiments.

As convergence criterion of the VBA-SR we use m k -m k-1 / m k-1 < 10 -5 , where m k and m k-1

represents the estimate of the HR image at kth and (k -1)th iteration, respectively. And for EGrad-SR and our approach, they stop when they achieve a HR image of a PSNR value very close (difference within 1‰) to that obtained by VBA-SR. . This result suggests that the introduction of the memory gradient subspace gives better directions and the proposed approach thus needs less iterations than the approach based on the gradient direction. Even though each iteration of the proposed approach takes more computation time than EGrad-SR due to its complexity, the proposed approach is still more efficient than EGrad-SR in most cases thanks to the decrease in the number of iterations. Concerning the VBA-SR, it takes less iterations than the other two approaches. However, each of its iteration takes much more time since it contains an inner loop due to the use of the conjugate gradient method to avoid the direct matrix inversion. As a result, VBA-SR is less efficient than the proposed approach in terms of computation time.

In order to compare the visual quality of reconstructed images, we show in Fig. 1 and Fig. 2 1 and Fig. 2, we can see that all the approaches increase the image resolution. Even in highly noisy cases where SNR = 5 dB, the noise present in LR images is effectively reduced meanwhile the image edges are not over-smoothed. Moreover, we can see that all the approaches give very similar HR images in each case, which is coherent with the similar PSNR values achieved by these approaches. Furthermore, as stated above, all the HR images were obtained with a same value of θ. We can see here that our approaches work well with this value in all the tested cases. This approach is applied to a linear inverse problem where a TV image prior with an approximate partition function and Jeffrey's hyperpriors are used, which results in a fully automatic algorithm. We have shown on a super-resolution problem that the proposed algorithm is much more fast than state of art approaches.

May 12, 2014 DRAFT

The reason is that we have integrated the memory gradient subspace optimization method which allows more flexibility in the choice of directions.

  which span the subspace and the vector s k = [s k 1 , ..., s k I ] T encloses the step-sizes along each direction. The subspace optimization method offers more flexibility in the choice of the descent direction d k by taking linear combinations of directions in D k .

1 , s subopt 2

 12 ), (54) b. Determine the suboptimal step-sizes (s subopt

  one of the LR images, the reconstructed images obtained by VBA-SR, EGrad-SR and our proposed approach for Cameraman and Jetplane LR images of SNR = 5, 15, 25, 35 dB. Comparing the LR images (shown in the top row) and the reconstructed HR images (given in the second, third and the bottom row) in Fig.

Fig. 1 :Fig. 2 :

 12 Fig. 1: One of the LR images (top row), HR images obtained by VBA-SR (the second row), EGrad-SR (the third row), Proposed approach (bottom row) for the Cameraman in the cases where SNR = (a) 5 dB; (b) 15 dB; (c) 25 dB (d) 35 dB. All images are presented in the same range of grayscale.

TABLE I :

 I PERFORMANCE COMPARISON OF VBA-SR[START_REF] Babacan | Variational Bayesian super resolution[END_REF], EGRAD-SR[START_REF] Fraysse | A measure-theoretic variational Bayesian algorithm for large dimensional problems[END_REF] AND OUR PROPOSED APPROACH IN TERMS OF NUMBER OF ITERATIONS/CPU TIME (IN SECONDS). TM) i7-3770 CPU (3.40 GHz) with 8.0 GB RAM. In TableI, we use bold numbers to indicate the best results, i.e. the shortest computation time for each data. Comparing the computation time, we can see that the proposed approach is much more efficient than the VBA-SR for all the datasets and more efficient than EGrad-SR in most cases, especially in lower noise ones such as cases where SNR = 25, 35 and 45 dB. For instance, in the case where SNR = 5 dB for the Cameraman image, the proposed approach takes 4.9 seconds which is only 9% of 56.9 seconds used by VBA-SR and slightly smaller than 5.0 seconds taken by EGrad-SR, in the case where SNR = 25 dB for the Lena image, the proposed approach took 5.4 seconds which is nearly 32% of the time needed by VBA-SR (16.8 seconds) and 71% of the time taken by EGrad-SR (7.6 seconds), in the case where SNR = 45 dB for the Jetplane image, the proposed approach took 7.9 seconds which is only 20% of the time taken by VBA-SR (40 seconds) and 43% of the time used by EGrad-SR (18.2 seconds). On average, the proposed approach is approximately 4 times faster than the VBA-SR and approximately 1.7 times as fast as EGrad-SR. Comparing the number of iterations, the proposed approach generally takes less iterations than the EGrad-SR (only one exception in the case where SNR = 15 dB of Cameraman image)

	Core(
	Data	PSNR	VBA-SR	EGrad-SR Proposed
	Testpat

May 12, 2014 DRAFT

h ∈ L 1 (µ) ⇔ R J |h(w)|µ(dw) < ∞