Kamélia Daudel
email: kamelia.daudel@telecom-paris.fr

Randal Douc
email: randal.douc@telecom-sudparis.eu

Franc ¸ois

THE F -DIVERGENCE EXPECTATION ITERATION SCHEME

Keywords: MSC 2010 subject classifications: Primary 60K35, 60K35; secondary 60K35 sample, L A

This paper introduces the f -EI(φ) algorithm, a novel iterative algorithm which operates on measures and performs f -divergence minimisation in a Bayesian framework. We prove that for a rich family of values of (f, φ) this algorithm leads at each step to a systematic decrease in the f -divergence and show that we achieve an optimum. In the particular case where we consider a weighted sum of Dirac measures and the α-divergence, we obtain that the calculations involved in the f -EI(φ) algorithm simplify to gradient-based computations. Empirical results support the claim that the f -EI(φ) algorithm serves as a powerful tool to assist Variational methods.

1. Introduction. Bayesian statistics for complex models often induce intractable and hard-to-compute posterior densities which need to be approximated. Variational methods such as Variational Inference (VI) [START_REF] Michael | An introduction to variational methods for graphical models[END_REF][START_REF] James | Variational algorithms for approximate bayesian inference[END_REF] and Expectation Propagation (EP) [START_REF] Opper | Gaussian processes for classification: Mean-field algorithms[END_REF][START_REF] Thomas | Expectation propagation for approximate bayesian inference[END_REF] consider this objective purely as a (non-convex) optimisation problem. These approaches seek to approximate the posterior density by a variational density q θ , characterized by a set of variational parameters θ ∈ T, where T is the parameter space. In these methods θ is optimised such that it minimizes a certain divergence D between the posterior and the variational density, typically the Kullback-Leibler (KL) divergence [START_REF] Kullback | On information and sufficiency[END_REF].

Modern optimisation-based approximate inference methods improved in three major directions [START_REF] Blei | Variational Inference: A Review for Statisticians[END_REF][START_REF] Zhang | Advances in Variational Inference[END_REF]. Firstly, Variational methods used to be limited to conditionally conjugate exponential family models [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF]. Monte Carlo methods and Black-Box inference techniques such as [START_REF] Paisley | Variational Bayesian Inference with Stochastic Search[END_REF][START_REF] Ranganath | Black Box Variational Inference[END_REF] have since been deployed, rendering Variational methods applicable to a wide range of models. These methods use gradients which are computed through automatic differentiation tools and climb the Monte Carlo approximated Variational Bound to the log-likelihood.

Secondly, classical Variational methods focused on the KL divergence as an objective function. However, in the VI case, the KL is known to often underestimate the variance and may miss important local modes of the true posterior [START_REF] Hoffman | Stochastic variational inference[END_REF][START_REF] Blei | Variational Inference: A Review for Statisticians[END_REF]. As for the EP algorithm, which performs local minimization, it is not guaranteed to converge and does not provide an easy estimate of the marginal likelihood [START_REF] Thomas | A Family of Algorithms for Approximate Bayesian Inference[END_REF]. Modern research consequently turned to more flexible families of divergence in order to reach better accuracy, the α-divergence [START_REF] Zhu | Bayesian invariant measurements of generalization[END_REF][START_REF] Zhu | Information geometric measurements of generalisation[END_REF] and Renyi's α-divergence [START_REF] Rényi | On measures of entropy and information[END_REF][START_REF] Van Erven | Rényi Divergence and Kullback-Leibler Divergence[END_REF] being such examples [START_REF] Minka | Divergence measures and message passing[END_REF][START_REF] Minka | Power ep[END_REF][START_REF] Miguel Hernández-Lobato | Black-box α-divergence Minimization[END_REF][START_REF] Li | Rényi Divergence Variational Inference[END_REF][START_REF] Bousso | Variational inference via χ upper bound minimization[END_REF].

Thirdly, scalable methods such as Stochastic Variational Inference [START_REF] Hoffman | Stochastic variational inference[END_REF] or Stochastic Expectation Propagation [START_REF] Li | Stochastic expectation propagation[END_REF][START_REF] Dehaene | Expectation Propagation in the large-data limit[END_REF] have been developed to enable large scale learning. These methods rely on stochastic optimisation techniques [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF][START_REF] Robbins | A stochastic approximation method[END_REF] and have been applied to complex probabilistic models, e.g Latent Dirichlet Allocation [START_REF] Blei | Latent dirichlet allocation[END_REF].

Framework. In this paper, we contribute in these three directions. The divergence we choose to work with is the f -divergence [START_REF] Morimoto | Eine informationstheoretische ungleichung und ihre anwendung auf den beweis der ergodizitat von markoffschen ketten[END_REF][START_REF] Morimoto | Markov processes and the h-theorem[END_REF], as it is a general family of divergences that encompasses the Kullback-Leibler, the reverse Kullback-Leibler and the α-divergence. Furthermore, we offer to change the space on which the minimization occurs. While the common approach consists in minimizing over the set of densities {y → q θ (y) : θ ∈ T} we consider instead a minimization over y → T µ(dθ)q θ (y) : µ ∈ M where M is a convenient subset of M 1 (T), the set of probability measures on T (and in this case, we equip T with a σ-field denoted by T ). In doing so, we extend the minimizing set to a larger space since a parameter θ can be identified with its associated Dirac measure δ θ . The paper is then organised as follows:

• In Section 2, we briefly review basic concepts around the f -divergence. We formulate explicitly the general form of the optimisation problem we consider, which includes the particular case of f -divergence posterior density approximation.

• In Section 3, we provide a new iterative scheme which we call the f -EI(φ) algorithm that performs an update of measures. We establish sufficient conditions on (f, φ) for this algorithm to lead at each step to a systematic decrease in the f -divergence and we obtain its convergence to an optimum. As the exact f -EI(φ) algorithm involves an integral which might be intractable, we define its approximate version and show its convergence to the exact algorithm.

• In Section 4, we apply the f -EI(φ) algorithm to f -divergence density approximation. In the particular case of the α-divergence and when µ is chosen as a weighted sum of Dirac measures, we obtain that the computations involved in the f -EI(φ) algorithm mostly rely on gradient-based calculations.

• Finally, Section 5 is devoted to numerical experiments, where we explore the impact of the hyperparameter φ as well as the state space dimension d and the parameter space dimension J on the convergence of the f -EI(φ) algorithm.

2. Formulation of the optimisation problem. Let (Y, Y, ν) be a measured space, where ν is a σ-finite measure on (Y, Y) and let f be a convex function over (0, ∞) that satisfies f (1) = 0. We start by defining the f -divergence between two probability measures P 1 and P 2 .

Definition 1 (f -divergence). Let P 1 and P 2 be two probability measures on (Y, Y) that are absolutely continuous with respect to ν i.e. P 1 ν, P 2 ν. Let us denote by p 1 = dP 1 dν and p 2 = dP 2 dν the Radon-Nikodym derivatives of P 1 and P 2 with respect to ν. The f -divergence between P 1 and P 2 is defined as follows :

D f (P 1 ||P 2 ) = Y f p 1 (y) p 2 (y) p 2 (y)ν(dy) . (1)
In [START_REF] Michael | An introduction to variational methods for graphical models[END_REF], we adopt the conventional notation 0f ( 0 0 ) = 0 and 0f

( a 0 ) = lim t↓0 tf ( a t ) = a lim u→∞ f (u)
u for all a > 0. We also extend the definition of f at 0 by setting f (0

) = lim t↓0 f (t) ∈ (-∞, ∞].
We now recall without proof a few results around the f -divergence and we refer to [START_REF] Cichocki | Families of alpha-beta-and gamma-divergences: Flexible and robust measures of similarities[END_REF][START_REF] Del | On asymptotic properties of informationtheoretic divergences[END_REF] and [START_REF] Sason | On f-divergences: Integral representations, local behavior, and inequalities[END_REF] for more details on the f -divergence family.

Proposition 2. The f -divergence is always non-negative and for f strictly convex, it is equal to zero if and only if P 1 = P 2 . Furthermore, it is jointly convex in P 1 and P 2 and

D f (P 1 ||P 2 ) = D f (P 2 ||P 1 ), where f (u) = uf ( 1 u ) is called the conjugate function.
Special cases include the Kullback-Leibler (KL) divergence D KL , the reverse Kullback-Leibler divergence D rKL and the α-divergence D (α)

A , where α ∈ R \ {0, 1} (see Table 1). The Hellinger distance [START_REF] Hellinger | Neue begründung der theorie quadratischer formen von unendlichvielen veränderlichen[END_REF][START_REF] Lindsay | Efficiency versus robustness: The case for minimum hellinger distance and related methods[END_REF] and the χ 2divergence [START_REF] Bousso | Variational inference via χ upper bound minimization[END_REF] correspond respectively to order α = 0.5 and α = 2 of the α-divergence. Moreover, the definition of the α-divergence can be extended to α ∈ {0, 1} by continuity and we recover the Kullback-Leibler when α → 1

f (u) Corresponding divergence u log(u) DKL(P1||P2) = Y log p 1 (y) p 2 (y) p1(y)ν(dy) -log(u) DrKL(P1||P2) = Y -log p 1 (y) p 2 (y) p2(y)ν(dy) 1 α(α -1) [u α -1] D (α) A (P1||P2) = 1 α(α-1) Y p 1 (y) p 2 (y) α p2(y)ν(dy) -1
and the reverse Kullback-Leibler when α → 0. Now consider a measurable space (T, T ). Let p be a measurable positive function on (Y, Y) and Q : (θ, A) → A q(θ, y)ν(dy) be a Markov transition kernel on T × Y with kernel density q defined on T × Y. Moreover, for all y ∈ Y, we denote µq(y) = T µ(dθ)q(θ, y) and we define, for all µ ∈ M 1 (T),

Ψ (f ) (µ) = Y f µq(y) p(y) p(y)ν(dy) . (2)
Note that p, q and ν appear as well in Ψ (f ) (µ) i.e Ψ (f ) (µ) = Ψ (f ) (µ; p, q, ν), but we drop them for notational ease and when no ambiguity occurs. Notice also that we replaced q θ (y) by q(θ, y) to comply with usual kernel notation.

We consider in what follows the general optimisation problem

arginf µ∈M Ψ (f ) (µ) , (3)
where p is a measurable positive function on (Y, Y).

This framework includes the particular case of f -divergence posterior density approximation. Indeed, let us consider the posterior density of the latent variables y given the data D: At this stage, a first remark is that the convexity of Ψ (f ) is straightforward from the convexity of f . Therefore, a simple yet powerful consequence of enlarging the variational family is that the optimisation problem now involves the convex mapping

µ → Ψ (f ) (µ) = Y f µq(y) p(y) p(y)ν(dy) ,
whereas the initial optimisation problem was associated to the mapping θ → Y f q θ (y) p(y) p(y)ν(dy), which is not necessarily convex. We now move on to Section 3, where we describe the f -EI(φ) algorithm and state our main theoretical results.

The f -Expectation Iteration algorithm f -EI(φ) .

3.1. An iterative algorithm for optimising Ψ (f ) . The following set of assumptions will be in force throughout the paper. We first gather the assumptions on q, p and ν.

(A1) The density kernel q on T × Y, the function p on Y and the σ-finite measure ν on (Y, Y) satisfy, for all (θ, y) ∈ T × Y, q(θ, y) > 0, p(y) > 0 and Y p(y)ν(dy) < ∞.

The next assumption mostly retains the assumptions on the function f used to define Ψ (f ) in ( 2).

(A2) The function f : (0, ∞) → R is monotonous, strictly convex and continuously differentiable, and f (1) = 0.

Under (A1) and (A2), we immediately obtain a lower bound on Ψ (f ) .

Lemma 3. Suppose that (A1) and (A2) hold. Then, for all µ ∈ M 1 (T), we have

Ψ (f ) (µ) f Y p(y)ν(dy) > -∞ ,
where f is defined in Proposition 2.

Proof. Since f (u) = uf (1/u), we have

Ψ (f ) (µ) = Y f p(y) µq(y) µq(y)ν(dy) .
Recalling that f and hence f , is convex on R >0 , Jensen's inequality applied to f yields

Ψ (f ) (µ) f Y p(y)ν(dy) > -∞.
Remark 4. Assumption (A1) can be extended by discarding the assumption that p(y) is positive for all y ∈ Y. As it complicates the expression of the constant appearing in the bound without increasing dramatically the degree of generality of the results, we chose to maintain this assumption for the sake of simplicity.

Thus, if there exists a sequence of probability measures {µ n : n ∈ N} on (T, T ) such that Ψ (f ) (µ 0 ) < ∞ and Ψ (f ) (µ n ) is non-increasing with n, Lemma 3 guarantees that this sequence converges to a limit in R. We now focus on constructing such a sequence {µ n : n ∈ N}.

For this purpose, let φ ∈ R := R \ {0}. The one-step transition of the f -Expectation Iteration algorithm f -EI(φ) can be formally described as an expectation step and an iteration step:

Algorithm 1: Exact f -EI(φ) one-step transition 1. Expectation step : bµ(θ) = Y q(θ, y)f µq (y) p(y) ν(dy)
2. Iteration step :

I φ (µ)(dθ) = µ(dθ) • |bµ(θ)| φ µ(|bµ| φ )
Given any initial measure µ ∈ M 1 (T) such that Ψ (f ) (µ) < ∞, the iterative sequence of probability measures (µ n ) n∈N is then defined by setting (4)

µ 0 = µ, µ n+1 = I φ (µ n ) , n ∈ N .
Note that under (A1) and (A2), b µ is well-defined (since f is of constant sign) and |b µ | ∈ (0, ∞] for all µ ∈ M 1 (T). Given φ ∈ R , the iteration µ → I φ (µ) is thus well-defined if moreover we have

0 < µ(|b µ | φ ) < ∞ . (5)
In the following part, we investigate some core properties of the aforementioned sequence of probability measures (µ n ) n∈N . For all µ ∈ M 1 (T) such that Ψ (f ) (µ) < ∞, we identify couples (f, φ) such that (i) f satisfies (A2), (ii) the sequence (µ n ) n∈N is well-defined and (iii) the f -EI(φ) algorithm diminishes Ψ (f ) (µ n ) at each iteration.

Monotonicity.

To establish that the f -EI(φ) algorithm diminishes Ψ (f ) (µ n ) at each iteration, we first derive a general upper-bound for the difference Ψ (f ) (ζ) -Ψ (f ) (µ). Here, (ζ, µ) is a couple of probability measures where ζ is dominated by µ, which we denote by ζ µ.

This first result relies on the existence of an exponent satisfying condition (A3) below, which will later on be used to specify a range of values for φ satisfying (5) for any µ ∈ M 1 (T).

(A3) We have ∈ R\[0, 1] and the function f : u → f (u 1/ ) is non-decreasing and concave on R >0 .

Proposition 5. Assume (A1), (A2) and (A3). Then, for all µ, ζ ∈ M 1 (T) such that µ(|b µ |) < ∞ and ζ µ,

(6) Ψ (f ) (ζ) Ψ (f ) (µ) + | | -1 {µ (|b µ |g ) -µ(|b µ |)} .
where g is the density of ζ wrt µ, i.e. ζ(dθ) = µ(dθ)g(θ). Moreover, equality holds in (6) if and only if ζ = µ.

Proof. Under (A2) f is of constant sign and under (A3), if > 1, f is non-decreasing and if < 0, f is non-increasing. This implies that for all u > 0 and all = 0, -1 f (u) = | | -1 |f (u)|, which will be used later in the proof.

Write by definition of f in (A3) and ζ, where the last inequality follows from Jensen's inequality applied to the convex function u → u (since ∈ R \ [0, 1]) and the fact that f is nondecreasing. Now set

Ψ (f ) (ζ) = Y f ζq (
u y = T µ(dθ) q(θ, y) µq(y) g(θ)µq(y) p(y) v y = µq(y) p(y)
and note that

u y -v y = µq(y) p(y) T µ(dθ) q(θ, y) µq(y) g (θ) -1 (8) Since f is concave, f (u y ) f (v y ) + f (v y )(u y -v y ).
Then, combining with [START_REF] Zhang | Advances in Variational Inference[END_REF], we get

Ψ (f ) (ζ) Y f (u y )p(y)ν(dy) (9) Y f (v y )p(y)ν(dy) + Y f (v y )(u y -v y )p(y)ν(dy)
Note that the first term of the rhs can be written as

Y f (v y )p(y)ν(dy) = Y f µq(y) p(y) p(y)ν(dy) = Ψ (f ) (µ) (10) Using now f (v y ) = -1 v 1/ -1 y f (v 1/
y ) and ( 8), the second term of the rhs of ( 9) may be expressed as

Y f (v y )(u y -v y )p(y)ν(dy) = -1 Y µq(y) p(y) 1- f µq(y) p(y) µq(y) p(y) T µ(dθ) q(θ, y) µq(y) g (θ) -1 p(y)ν(dy) = -1 T µ(dθ) Y q(θ, y)f µq(y) p(y) ν(dy) g (θ) --1 Y µq(y)f µq(y) p(y) ν(dy) = | | -1 {µ (|b µ |g ) -µ(|b µ |)} .
Combining this equality with ( 9) and ( 10) finishes the proof of the inequality.

If the equality holds in ( 6), then the equality in Jensen's inequality [START_REF] Zhang | Advances in Variational Inference[END_REF] shows that g is constant µ-a.e. so that ζ = µ, and the proof is completed.

We now plan on setting ζ = I φ (µ) in Proposition 5 and obtain that one iteration of the f -EI(φ) algorithm yields Ψ (f ) • I φ (µ) Ψ (f ) (µ). For this purpose and based on the upper bound obtained in Proposition 5, we strengthen the condition (5) as follows to take into account the exponent

0 < µ(|b µ | φ ) < ∞ and µ(|b µ |g ) µ(|b µ |) with g = |bµ| φ µ(|bµ| φ ) (11)
This leads to our first main theorem.

Theorem 1. Assume (A1), (A2) and (A3). Let µ ∈ M 1 (T) be such that µ(|b µ |) < ∞ and let φ ∈ R satisfy [START_REF] Hoffman | Stochastic variational inference[END_REF]. Then, the two following assertions hold.

(i) We have Ψ (f ) • I φ (µ) Ψ (f ) (µ). (ii) We have Ψ (f ) • I φ (µ) = Ψ (f ) (µ) if and only if µ = I φ (µ).
Proof. We apply Proposition 5 with

ζ = I φ (µ) so that ζ(dθ) = µ(dθ)g(θ) with g = |b µ | φ /µ(|b µ | φ ). Then, (12) Ψ (f ) • I φ (µ) Ψ (f ) (µ) + | | -1 {µ (|b µ |g ) -µ(|b µ |)} Ψ (f ) (µ)
where the last inequality follows from condition [START_REF] Hoffman | Stochastic variational inference[END_REF].

Let us now show (ii). The if part is obvious. As for the only if part,

Ψ (f ) • I φ (µ) = Ψ (f ) (µ) combined with (12) yields Ψ (f ) • I φ (µ) = Ψ (f ) (µ) + | | -1 {µ (|b µ |g ) -µ(|b µ |)} ,
which is the case of equality in Proposition 5. Therefore, I φ (µ) = µ.

We are now able to derive our second main theorem.

Theorem 2. Assume that p and q are as in (A1). Let (f, φ) belong to any of the following cases.

(i) Reverse Kullback-Leibler: f (u) = -log(u), and φ ∈ (0, 1].

(ii) α-Divergence: f (u) = 1 α(α-1) (u α -1), (a) α ∈ (-∞, -1] and φ ∈ (0, -1/α]; (b) α ∈ (-1, 1) \ {0} and φ ∈ (0, 1]; (c) α ∈ (1, ∞) and φ ∈ (1/(1 -α), 0). Then (A2) holds. Moreover, let µ ∈ M 1 (T) be such that Ψ (f ) (µ) < ∞.
Then the sequence (µ n ) n∈N defined by (4) is well-defined and the sequence

(Ψ (f ) (µ n )) n∈N is non-increasing.
The proof of this theorem requires intermediate results, which are derived in Appendix A.1 alongside with the proof of Theorem 2. Remark 6. In the proof of Theorem 2, we also prove along the way that either

Ψ (f ) (µ n+1 ) < Ψ (f ) (µ n ) for all n 0, or that there exists n 0 < ∞ such that µ n = µ n 0 for all n n 0 i.e µ n 0 = I φ (µ n 0 ).
The results we obtained at this point are summarized in Table 2. We now study the limiting behavior of the f -EI(φ) algorithm for the iterative sequence of probability measure (µ n ) n∈N .

Table 2

Allowed ranges for φ in the f -EI(φ) algorithm per divergence

Divergence considered

Corresponding range

Reverse KL f (u) = -log(u) φ ∈ (0, 1] α-divergence f (u) = 1 α(α-1) (u α -1) α ∈ (-∞, -1] φ ∈ (0, -1/α] α ∈ (-1, 1) \ {0} φ ∈ (0, 1] α ∈ (1, ∞) φ ∈ (1/(1 -α), 0) 3.3.
Limiting behavior of the Exact f -EI(φ) algorithm. Let µ ∈ M 1 (T) and let us consider the iterative sequence of probability measures (µ n ) n∈N defined by [START_REF] Thomas | Expectation propagation for approximate bayesian inference[END_REF]. We seek to identify the possible weak limits of (µ n ) n∈N . To do so, we introduce the following additive set of assumptions (A4) (i) T is a compact metric space and T is the associated Borel σ-field;

(ii) for all y ∈ Y, θ → q(θ, y) is continuous;

(iii) we have

Y sup θ∈T f q(θ, y) p(y) p(y)ν(dy) < ∞;
(iv) we have

Y sup θ∈T q(θ, y) × sup θ ∈T f q(θ , y) p(y) ν(dy) < ∞.
Notice that conditions (A4)-(iii) and (A4)-(iv) act as a uniform majoration of Ψ (f ) (µ) and b µ (θ) with respect to µ and θ. In addition, these two conditions are rather weak under (A4)-(i), since we consider in each a supremum taken over a compact set (and T will always be chosen as such in practice).

We define M 1,µ (T) as the set of probability measures dominated by µ. We have the following theorem, which states that the possible weak limits of (µ n ) n∈N corresponds to the global infimum of Ψ (f ) . Theorem 3. Assume (A1), (A2) and (A4). Then, for all ζ ∈ M 1 (T) any φ ∈ R satisfies (5) and

Ψ (f ) (ζ) < ∞.
Further assume that there exists µ, μ ∈ M 1 (T) such that the (well-defined) sequence (µ n ) n∈N defined by (4) weakly converges to μ as n → ∞. Then the following assertions hold

(i) μ is a fixed point of I φ , (ii) Ψ (f ) (μ) = inf ζ∈M 1,µ (T) Ψ (f ) (ζ),
in any of the following cases:

(a) f non-increasing and φ > 0, (b) f non-decreasing and φ < 0.

The proof of Theorem 3 can be found in Appendix A.2. Note that any couple (f, φ) described in Table 2 falls into one of the two categories (a) f non-increasing and φ > 0 or (b) f non-decreasing and φ < 0.

Remark 7. In the particular case of the α-divergence, for which

f (u) = 1 α(α-1) [u α -1], Condition (A4)-(iv) can be rewritten as Y sup θ∈T q(θ, y) × sup θ ∈T q(θ , y) p(y) α-1 ν(dy) < ∞ ,
which in turns, implies that Condition (A4)-(iii) is satisfied.

Approximate f -EI(φ) .

As Algorithm 1 typically involves an intractable integral in the Expectation step, we now turn to a practical version of this algorithm.

Algorithm 2: Approximate f -EI(φ) one-step transition 1. Sampling step : Draw independently Y1, ..., YK ∼ µq 2. Expectation step : bµ,K (θ) = 1 K K k=1 q(θ, Y k ) µq(Y k ) f µq(Y k ) p(Y k ) 3. Iteration step : I φ K (µ)(dθ) = µ(dθ) • |bµ,K (θ)| φ µ(|bµ,K | φ )
Algorithm 2 uses µq as a sampler instead of q(θ, •). Indeed, as our algorithm optimises over µ, sampling with respect to µq gives preference to the interesting regions of the parameter space. Furthermore, picking a sampler that is independent from θ is less costly from a computational point of view.

In the rest of this section, we consider i.i.d random variables Y 1 , Y 2 , ... with common density µq w.r.t ν, defined on the same probability space (Ω, F, P) and we denote by E the associated expectation operator.

Proposition 8. Assume (A1) and (A2). Let µ ∈ M 1 (T), φ ∈ R be such that µ(|b µ |) ∨ µ(|b µ | φ ) < ∞ and (13) T µ(dθ)E q(θ, Y 1 ) µq(Y 1 ) f µq(Y 1 ) p(Y 1 ) φ < ∞ .
Then, lim

K→∞ I φ K (µ) -I φ (µ) T V = 0 , P -a.s.
Proof. By the triangular inequality, for all K ∈ N , for all θ ∈ T,

|b µ,K (θ)| φ µ(|b µ,K | φ ) - |b µ (θ)| φ µ(|b µ | φ ) |b µ,K (θ)| φ µ(|b µ,K | φ ) 1 - µ(|b µ,K | φ ) µ(|b µ | φ ) + ||b µ,K (θ)| φ -|b µ (θ)| φ | µ(|b µ | φ )
Thus,

I φ K (µ) -I φ (µ) T V = µ |b µ,K | φ µ(|b µ,K | φ ) - |b µ | φ µ(|b µ | φ ) 1 - µ(|b µ,K | φ ) µ(|b µ | φ ) + µ(||b µ,K | φ -|b µ | φ |) µ(|b µ | φ )
For the first term of the rhs, Lemma 21 yields

lim K→∞ 1 - µ(|b µ,K | φ ) µ(|b µ | φ ) = 0 (14)
As for the second term of the rhs, first note that for all K ∈ N , for all θ ∈ T Lemma 21 and (16)

0 ||b µ,K (θ)| φ -|b µ (θ)| φ | |b µ,K (θ)| φ + |b µ (θ)| φ , ( 15
)
and since µ(|b µ | φ ) < ∞ the LLN for µ-almost all θ ∈ T yields lim K→∞ |b µ,K (θ)| φ = |b µ (θ)| φ . (16) Furthermore, since µ(|b µ | φ ) < ∞,
imply lim K→∞ µ |b µ,K | φ + |b µ | φ = µ lim K→∞ |b µ,K | φ + |b µ | φ < ∞
Combining with ( 15) and ( 16), we apply Lemma 22 and obtain

lim K→∞ µ(||b µ,K | φ -|b µ | φ |) µ(|b µ | φ ) = 0
which, along with [START_REF] Zhu | Information geometric measurements of generalisation[END_REF], finishes the proof.

Note that since Algorithm 2 is to be repeated multiple times until the convergence is reached, Ψ (f ) (µ) needs to be approximated as well, as it measures the convergence. In particular, the expression of Ψ (f ) (µ) given in (2) legitimates the choice of µq as the sampler, as the closest available approximation of p. Let us denote by Ψ (f )

K (µ) the corresponding unbiased approximation of Ψ (f ) (µ) Ψ (f ) K (µ) = 1 K K k=1 f µq(Y k ) p(Y k ) p(Y k ) µq(Y k ) .
Lemma 9. Assume (A1) and that (f, φ) belongs to any of the cases stated in Table 2

. Let µ ∈ M 1 (T) be such that Ψ (f ) (µ) < ∞. Then for all K ∈ N , Ψ (f ) K • I φ K (µ) Ψ (f ) K (µ) .
Furthermore,

lim K→∞ Ψ (f ) K (µ) = Ψ (f ) (µ) , P -a.s.
Proof. The first point is a straightforward adaptation of the proof of Theorem 2.

As for the second point, we know from Lemma 3 that

E[|f ( µq(Y 1 ) p(Y 1 ) )| p(Y 1 ) µq(Y 1 )
] is finite if and only if Ψ (f ) (µ) is finite. We can thus apply the LLN which yields the desired result.

In the next section, we investigate how our algorithm can be applied to density approximation for the divergences identified in Theorem 2.

f -EI(φ) applied to density approximation.

4.1. Reformulation of the optimisation problem. Let p be a probability density function on (Y, Y) and assume that we only have access to an unnormalized version p * of the density p, that is for all y ∈ Y,

p(y) = p * (y) Z , (17)
where Z := Y p * (y)ν(dy) is called the normalizing constant or partition function.

Let us denote by P the probability measure on (Y, Y) with density p with respect to ν and let us recall that for all µ ∈ M 1 (T), µQ corresponds to the probability measure on (Y, Y) with density µq with respect to ν. We then have the following lemma, whose proof can be found in Appendix C. Lemma 10. Assume (A1). Then, for both the reverse Kullback-Leibler and the α-divergence, optimising the objective D f (µQ|| P) (with respect to µ) is equivalent to optimising the objective Ψ (f ) (µ; p) with p = p * .

If we now rewrite Lemma 3 for p = p * in the particular case of the reverse Kullback-Leibler and of the α-divergence, we obtain, as f (u) = u log(u) and

f (u) = 1 α(α-1) [u 1-α -u] respectively, Y -log µq(y) p * (y) p * (y)ν(dy) Z log(Z) , (18)
Y 1 α(α -1) µq(y) p * (y) α p * (y)ν(dy) 1 α(α -1) Z 1-α . ( 19
)
Here, the normalizing constant Z only appears in the r.h.s of ( 18) and ( 19). An interesting aspect is that optimising Ψ (f ) (µ; p) with p = p * is equivalent to optimising the bound [START_REF] Minka | Power ep[END_REF] for the reverse Kullback-Leibler and the bound [START_REF] Miguel Hernández-Lobato | Black-box α-divergence Minimization[END_REF] for the α-divergence, where the optimisation does not involve the (unknown) normalizing constant Z anymore.

As it turns out, ( 18) is of little help to provide an explicit bound on the normalizing constant Z, however [START_REF] Miguel Hernández-Lobato | Black-box α-divergence Minimization[END_REF] can be used to provide either an upper or lower bound on Z, depending on the sign of α. To see this, let α ∈ R \ {0, 1}. Then, for any measurable positive function q on (Y, Y), set

ξ (α) (q) := Y q(y) p * (y) α p * (y)ν(dy) 1 1-α . ( 20
)
We call α-bound the function q → ξ (α) (q). The next lemma is a straightforward consequence of [START_REF] Miguel Hernández-Lobato | Black-box α-divergence Minimization[END_REF].

Lemma 11. Assume (A1). Let µ ∈ M 1 (T). Then, for all α + ∈ (0, 1) ∪ (1, +∞) and all α -< 0, we have

ξ (α + ) (µq) Z ξ (α -) (µq) .
Of course, the two previous lemmas go beyond the case q = µq. Notably, Lemma 11 holds for any probability density function q on (Y, Y) such that supp(p * (y)) ⊆ supp(q(y)), as previously established in [START_REF] Minka | Divergence measures and message passing[END_REF]Theorem 2]. Furthermore, α → ξ (α) (q) is continuous on α : ξ (α) (q) < +∞ , which is a straightforward consequence of [20, Theorem 1] using Renyi's α-divergence. In conclusion, we acquire additional information by using the α-divergence instead of the reverse Kullback-Leibler, since the α-bound can be used as a bound for the normalizing constant Z. Then, under the assumptions of Theorem 2, we can apply Algorithm 1 iteratively and benefit from the fact that (i) the α-bound acts as a measure of convergence and (ii) the monotonicity property can be observed through (ξ (α) (µ n q)) n∈N , as for all µ ∈ M 1 (T),

Ψ (f ) (µ; p) = ξ (α) (µq) 1-α - Z α(α -1)
with p = p * .

We now focus on rewriting Algorithm 2, which is the algorithm that we use in practice, in the special case of the α-divergence and with p = p * . 4.2. Approximate f -EI(φ) for the α-divergence. From here on, we only consider the particular case of the α-divergence with p = p * . Let (α, φ) be as in Table 2 and let µ ∈ M 1 (T). Using that f (u) = 1 α-1 u α-1 , we obtain Algorithm 3.

(θ) = 1 K(1 -α) K k=1 q(θ, Y k )µq(Y k ) α-2 p * (Y k ) 1-α 3. Iteration step : I φ K (µ)(dθ) = µ(dθ) • |bµ,K (θ)| φ µ(|bµ,K | φ )
Now recall that ξ (α) (µq) acts as a surrogate to Ψ (f ) (µ; p) which does not involve the normalizing constant Z anymore. If we are to apply Algorithm 3 repeatedly, we thus need to approximate this quantity in order to assess the convergence. Consequently, we define the Monte Carlo approximation ξ

(α) K (µq) of ξ (α) (µq) by ξ (α) K (µq) = 1 K K k=1 p * (Y k ) µq(Y k ) 1-α 1 1-α .
where Y 1 , ..., Y K are drawn independently from µq as in Algorithm 3. This estimator converges to ξ (α) (µq) as K → ∞ and while it is biased, the bias can be characterized using [START_REF] Li | Rényi Divergence Variational Inference[END_REF]Theorem 2].

Finally, Algorithm 3 requires us to know how to sample from µq. We address the case where µ corresponds to a weighted sum of Dirac measures. This case is of particular interest to us since, as we shall see, it provides a gradient-based update formula for the mixture weights. To this end, let J ∈ N and let us introduce the simplex of R J S J =    λ = (λ 1 , ..., λ J ) ∈ R J : ∀j ∈ {1, ..., J} , λ j 0 and

J j=1 λ j = 1    .
Let θ 1 , ..., θ J ∈ T be fixed. For all λ ∈ S J , we define µ λ ∈ M 1 (T) by

µ λ = J j=1 λ j δ θ j .
Then, for all λ ∈ S J , µ λ q(y) = J j=1 λ j q(θ j , y) corresponds to a mixture model. Let (µ n ) n∈N be defined by

µ 0 = µ λ , µ n+1 = I φ K (µ n ) , n ∈ N .
An immediate induction yields that for every n ∈ N, µ n can be expressed as µ n = J j=1 λ j,n δ θ j where λ n = (λ 1,n , . . . , λ J,n ) ∈ S J satisfies the initialisation λ 0 = λ and the update formula: for all n ∈ N and all j ∈ {1, . . . , J},

λ j,n+1 = λ j,n |b µn,K (θ j )| φ µ n (|b µn,K | φ ) , (21)
with Y 1,n , ..., Y K,n drawn independently from µ n q and b µn,K (θ

j ) = 1 K(1 -α) K k=1 q(θ j , Y k,n )µ n q(Y k,n ) α-2 p * (Y k,n ) 1-α .
We are thus able to derive Algorithm 4 below, where ( 21) is iterated until the convergence is reached and the α-bound is used at each step as a tractable measure of convergence and as a bound on the normalizing constant Z. Expectation step : Compute A λ = (aj) 1 j J where

aj = 1 K K k=1 q(θj, Y k )µ λ q(Y k ) α-2 p * (Y k ) 1-α (22)
and deduce B λ = (λja φ j ) 1 j J , b λ = J j=1 λja φ j and c λ = J j=1 λjaj.

Iteration step : Set ξ (α) K (µ λ q) ← c 1/(1-α) λ λ ← 1 b λ B λ end
In this particular framework, most of the computing effort at each step lies within the computation of the vector (b µn,K (θ j )) 1 j J , or equivalently the vector A λ = (a j ) 1 j J where the a j are defined as in [START_REF] Li | Stochastic expectation propagation[END_REF]. Interestingly, these computations are similar to the ones required in typical gradient-based variational methods involving the α-divergence or Renyi's α-divergence [START_REF] Miguel Hernández-Lobato | Black-box α-divergence Minimization[END_REF][START_REF] Li | Rényi Divergence Variational Inference[END_REF][START_REF] Bousso | Variational inference via χ upper bound minimization[END_REF]. Indeed, the objective function in these methods has a straightforward connection with the function q → L (α)

A (q) := Y 1 α(α -1) q(y) p * (y) α p * (y)ν(dy) ,
and in our case, the score gradient of L (α) A (µ λ q) is directly linked to the quantities that are approximated in the Mixture α-Approximate f -EI(φ) algorithm, since under the proper (differentiation) assumptions

∇ λ L (α) A (µ λ q) = (b µ λ (θ j )) 1 j J ,
where for all j ∈ {1, ..., J}, b µ λ (θ j ) = 1 α-1 Y q(θ j , y)

µ λ q(y) p * (y) α-1
ν(dy) .

Remark 13. Our (exact) mixture weights update rule states that for all n ∈ N, for all j ∈ {1, . . . , J},

λ j,n+1 = λ j,n |b µ λ n (θ j )| φ µ λ n (|b µ λ n | φ ) ,
In the particular case where (α, φ) = (0, 1) or (α, φ) = (-1, 1), notice that we recover respectively the weights update rules from the Population Monte Carlo algorithm for mixtures applied to reverse Kullback-Leibler minimisation [START_REF] Douc | Convergence of adaptive mixtures of importance sampling schemes[END_REF] and to Variance minimisation [START_REF] Douc | Minimum variance importance sampling via population monte carlo[END_REF].

We now move on to numerical experiments in the next section.

Numerical experiments.

To illustrate Algorithm 4, we first consider an example where p * corresponds to a mixture of two one-dimensional Gaussian densities multiplied by a positive constant Z, such that

p * (y) = Z × [γ 1 N (y; -s, 1) + γ 2 N (y; s, 1)] ,
where

γ 1 , γ 2 > 0, γ 1 + γ 2 = 1 and Z = 2.
This simple framework allows us to visualize and characterize the behavior of the algorithm with respect to φ. Looking back at [START_REF] Bousso | Variational inference via χ upper bound minimization[END_REF], we are prone to think that φ acts as a learning rate, as large values of |φ| amplify the impact of each |b µ,K (θ j )| in the update process while φ = 0 corresponds to no update at all. Let γ 1 = 0.8, γ 1 = 0.2, Q be a Gaussian kernel with variance 1, J = 2 with {θ 1 , θ 2 } = {-2, 2} and K = 5000. We perform 50 iterations of the α-Approximate f -EI(φ) algorithm. The results for α = -2 (with corresponding range (0, 0.5] for φ) can be seen on Figure 1 and several more examples are available in Appendix D.

As expected, we obtain faster convergence rates as φ gets bigger in absolute value. However, if |φ| is too large, the algorithm may oscillate around the true value of the parameters due to the discretization, as exemplified in Figure 1-(1). Thus, the hyperparameter φ behaves like a learning rate. Secondly, we are able to observe the monotonicity property from Theorem 2 in Figure 12, which plots the α-bound at each time step for various values of φ. Here, as we have picked α negative and as the conditions of support are met (supp(p * ) = supp(µ λ q)), Z is upper-bounded by the α-bound. Furthermore, since the true value of the parameters belongs to the optimisation set, the bound is attained as the optimisation is carried out. (1)

(2)

(3) (4) (5)
Towards an adaptive algorithm. Algorithm 4 leaves {θ 1 , ..., θ J } unchanged throughout the optimisation of the mixture weights (we call it an Exploitation Step). A natural idea is to combine this algorithm with an Exploration step that modifies the parameter set. We offer to derive the new parameter set from the old one using the following simple update rule: assume that the mixture weights have been optimised using Algorithm 4. We first resample among {θ 1 , ..., θ J } according to the optimised mixture weights. The obtained sample θ 1 , ..., θ J is then perturbed stochastically using a Gaussian transition kernel Q r with covariance matrix rI d (r > 0) and density q r , which gives us our new parameter set. The goal is then to iterate this Exploitation-Exploration procedure whilst using a decaying r (to ensure the convergence of the parameter set). As for the initial parameter set, it is generated randomly from an initial density q 0 , where we have in mind that q 0 should allow to explore the space extensively.

The complete algorithm is summed up in Algorithm 5. Note that we kept K and J fixed for convenience, but that they could be set according to a policy (K i ) i and (J i ) i and similarly we could consider a sequence of Markov transition kernels (Q i ) i . K (µ λ q). Exploration step : Draw independently J samples θ1, ..., θJ from µ λ qr i and set i = i + 1.

end

We now want to assess how Algorithm 5 performs in a higher dimensional setting. Notably, we aim at observing the impact of the dimension d as well as the parameter set dimension J. where u d is the d-dimensional vector whose coordinates are all equal to 1, I d is the identity matrix, s = 2 and Z = 2.

The first aspect that we need to consider is the choice of the divergence measure. The hyperparameter α allows to choose between mass-covering divergences which tends to cover all the modes (α 0) and mode-seeking divergences that are attracted to the mode with the largest probability mass (α 1), the case α ∈ (0, 1) corresponding to a mix of the two worlds. Depending on the learning task, the optimal α may differ and understanding how to select the value of α is still an area of ongoing research. In our case, since the targeted density is multimodal, we prefer having α < 1 and we set (α, φ) = (0.5, 1) in our experiments.

We set the initial sampler to be a centered normal distribution with covariance matrix 5I d , Q to be a Gaussian kernel with variance I d and (r i ) i which satisfies r 0 = 2.5 and r i = r 0 / √ i + 1 for all i ∈ N. Let K = 500, J ∈ {50, 100, 200, 400, 800}, d = {2, 4, 8, 16}. We run 10 iterations of Algorithm 5, with 25 inner iterations each time the Mixture α-Approximate f -EI(φ) algorithm is called. We replicate the experiment 100 times and calculate the average α-bound computed over the 100 replicates. The results can be seen on Figure 2.

Observe that the jumps in the α-bound correspond to an update of the parameter set. As expected, the optimisation becomes harder as the dimension grows. Yet, we are still able to climb to the normalizing constant Z up to the dimension 8 with J = 800, even though the Exploration step has not been optimised.

Conclusion and perspectives.

We introduced the f -EI(φ) algorithm, an iterative algorithm which operates on measures. We proved that for a rich family of values of (f, φ) this algorithm leads at each step to a systematic decrease in the f -divergence and obtained its convergence to an optimum. In the particular case of the α-divergence with µ set as a weighted sum of Dirac measures, we obtained that the mixture weights update rule mostly relied on gradient-based calculations. Empirical results confirmed that the hyperparameter φ acted as a learning rate for our algorithm. They also demonstrated that the f -EI(φ) algorithm serves as a powerful Exploitation step, which shall be combined with an appropriate Exploration step to form a fully adaptive algorithm.

To conclude, we state several directions to extend our work on both a theoretical and a practical level.

Learning rate. We maintained φ constant in the f -EI(φ) algorithm. Exploring variants of the algorithm with different decaying learning rate policies (φ n ) n and investigating convergence rates might result in more accurate and more stable results in practice.

Large scale learning. By noticing that the f -EI(φ) algorithm falls into the category of gradient-based algorithms when µ is chosen as a weighted sum of Diracs, we paved the way for large scale learning by using stochastic optimisation techniques, as deployed in [START_REF] Miguel Hernández-Lobato | Black-box α-divergence Minimization[END_REF] or more recently in [START_REF] Li | Rényi Divergence Variational Inference[END_REF][START_REF] Bousso | Variational inference via χ upper bound minimization[END_REF].

Exploration

Step. The f -EI(φ) algorithm allows us to extend the parameter set and to work with a population of particles {θ 1 , ..., θ J } instead of just one particle θ. In this regard, many more evolved methods could be envisioned as an Exploration step and combined with the f -EI(φ) algorithm.

Monte Carlo Approximation. One may want to resort to more advanced Monte Carlo methods in the estimation of b µn at each step. For example, we did not reuse any of the past samples so far in our calculations. Since we kept the allocation policy (K n ) n constant equal to K in this paper, another interesting aspect would be to investigate how different allocation policies affect the performances of the algorithm.

While Lemma 15 seems to advocate for g = dζ/dµ to be proportional to |b µ | 1/(1-) , notice that this choice of g might not be optimal to minimize Ψ (f ) (ζ) -Ψ (f ) (µ), as A(µ, ζ) also depends on g through ζ. In the next lemma, we thus propose another choice of the tuning parameter φ, which also satisfies [START_REF] Hoffman | Stochastic variational inference[END_REF]

for any µ ∈ M 1 (T) such that µ(|b µ |) < ∞. Lemma 16. Assume (A1), (A2) and (A3). Let µ ∈ M 1 (T) such that µ(|b µ |) < ∞. Assume in addition that | |
1, then the real number φ = -1/ satisfies [START_REF] Hoffman | Stochastic variational inference[END_REF].

Proof. Setting g ∝ |b µ | -1/ , we get µ(|b µ |g ) = µ(|b µ | 1-/ )[µ(|b µ | -1/ )] -= [µ(|b µ | -1/ )] - µ(|b µ |)
where the last inequality follows from Jensen's inequality applied to the convex function u → u -(since | | 1). Since µ(|b µ |) < ∞, the parameter φ = -1/ satisfies [START_REF] Hoffman | Stochastic variational inference[END_REF]. Furthermore using µ(|b µ |) < ∞ when φ < 0 and (A1) combined with (A2) for φ > 0, we obtain µ(|b µ | φ ) > 0, which concludes the proof.

Lemma 15 and Lemma 16 allow us to define a range of values for φ that decreases Ψ (f ) (µ n ) at each iteration step. Now, in order to prove Theorem 2, we need to check that the reverse Kullback-Leibler and the α-divergence are f -divergences with a function f that satisfies (A2) and (A3).

Proof of Theorem 2. The proof consists in verifying that we can apply Theorem 1, that is, for the two considered functions f , we must check (A2) and find a range of constants which satisfy (A3). We then use Lemma 15 or Lemma 16 to deduce that, for the provided constants φ, [START_REF] Hoffman | Stochastic variational inference[END_REF] holds for all µ n with n 0.

(i) Assumption (A2) readily holds and so does (A3) for all < 0, with f (u) = -log(u)/ . Moreover, by definition of b µn , we get for all n ∈ N,

µ n (|b µn |) = Y µ n q(y) p(y) µ n q(y) ν(dy) = Y p(y)ν(dy) < ∞ .
Combining with Lemma 15 and Lemma 16, [START_REF] Hoffman | Stochastic variational inference[END_REF] holds for all µ n with n 0 and for any φ ∈ (0, 1].

(ii) Again (A2) can be readily checked. Observing that for α / ∈ {0, 1},

f (u) = 1 α(α -1) u α/ -1 ,
we get that (A3) holds for

     α if α < 0, < 0 if α ∈ (0, 1), α if α > 1.
Lemmas 15 and 16 provide the corresponding ranges for φ in Cases (a), (b) and (c). To finish the proof, we now show by induction that for all n ∈ N, Ψ (f ) (µ n ) and µ n (|b µn |) are finite, so that Lemmas 15 and 16 can indeed be applied.

Since uf (u) = αf (u) + 1/(α -1), we have, for all n ∈ N,

µ n (|b µn |) = Y µ n q(y) p(y) f µ n q(y) p(y) p(y)ν(dy) (25) |α| Y f µ n q(y) p(y) p(y)ν(dy) + 1 |α -1|
Using Lemma 3, the rhs is finite if and only if Ψ (f ) (µ n ) is finite. We can now check that Ψ (f ) (µ n ) and µ n (|b µn |) are finite by induction on n:

• Start with n = 0. Then µ 0 = µ and by assumption, Ψ (f ) (µ) < ∞ and we deduce µ(|b µ |) < ∞ by [START_REF] Robbins | A stochastic approximation method[END_REF].

• Now, under the induction assumption and using Lemma 15 and Lemma 16, Theorem 1 can be applied and Ψ (f ) (µ n+1 ) Ψ (f ) (µ n ) < ∞. Using again [START_REF] Robbins | A stochastic approximation method[END_REF] with n replaced by n + 1, we get µ n+1 (|b µ n+1 |) < ∞.

A.2. Proof of Theorem 3. In the following, we use the notation µ n ⇒ μ for the weak convergence of measures in M 1 (T). We first derive four useful lemmas.

Lemma 17. Assume (A1), (A2) and (A4). Suppose that µ n ⇒ μ. Then the following assertions hold.

(i) For all y ∈ Y, µ n q(y) tends to μq(y) as n → ∞. Proof. We prove the assertions successively. Proof of (i). For all y ∈ Y, the function θ → q(θ, y) is continuous on a compact set, hence bounded. The weak convergence µ n ⇒ μ thus implies the pointwise convergence of µ n q to μq. Proof of (ii). For all θ ∈ T and ζ ∈ M 1 (T), we write

|a ζ (θ, y)| sup θ ∈T q(θ , y) × sup θ ∈T f q(θ , y) p(y) , (26)
which is integrable w.r.t ν(dy) by (A4)-(iv). The second part of (ii) is obtained similarly. Using (i) and that f is C 1 by (A2), we get that, for all (θ, y) ∈ T × Y, lim n→∞ q(θ, y) f µ n q(y) p(y) = q(θ, y) f μq(y) p(y) ,

i.e lim

m -:= Y inf θ ∈T q(θ , y) × inf θ ∈T f q(θ , y) p(y) ν(dy) , (27)
m + := Y sup θ ∈T q(θ , y) × sup θ ∈T f q(θ , y) p(y) ν(dy) .
We have that m + is finite by (A4)-(iv). Now recall that under (A2), f does not vanish on (0, ∞). Together with (A1), we thus have that for any y ∈ Y, the functions θ → q(θ, y) and θ → |f (q(θ, y)/p(y)) | are continuous and positive on the compact set T, from which we deduce that m -> 0.

Proof of (iv). . Now, for all y ∈ Y,

|a n (y)| 2 sup θ∈T q(θ, y) × sup θ ∈T f q(θ , y) p(y)
which is integrable w.r.t ν by (A4)-(iv). Moreover, by (i) and by continuity of f , we have lim n→∞ a n (y) = 0, and (29) follows by dominated convergence. Finally, combining ( 28), ( 29) and

µ n |b µn | φ h = µ n |b µn | φ h -|b μ| φ h + µ n |b μ| φ h ,
we obtain (iv), and the proof is concluded.

Lemma 18. Assume (A1) and (A2). Let ζ, ζ ∈ M 1 (T). Then,

T [ζ -ζ ](dθ)b ζ (θ) Ψ (f ) (ζ) -Ψ (f ) (ζ ) . (30)
Let μ, µ ∈ M 1 (T) and assume that there exists µ ∈ M 1,µ (T) such that Ψ (f ) (µ ) < Ψ (f ) (μ). Then, for f non-increasing (resp. non-decreasing), there exists δ > 1 (resp. δ < 1) and such that

µ (b μ < δ μ(b μ)) > 0 . (31) Proof. By definition of b ζ , T (ζ -ζ )(dθ)b ζ (θ) = T (ζ -ζ )(dθ) Y q(θ, y)f ζ q(y) p(y) ν(dy) = Y ζq(y) -ζ q(y) p(y) f ζ q(y) p(y) p(y)ν(dy) .
Assume that (µ n ) n∈N weakly converges to μ ∈ M 1 (T). First note that Lemma 17-(iii) implies that for any φ ∈ R the sequence (µ n ) n∈N is welldefined and μ satisfies [START_REF] Kullback | On information and sufficiency[END_REF].

We now prove Assertions (i) and (ii) successively. Proof of (i). For all ζ ∈ M 1 (T) and all y ∈ Y, set a ζ (y) = f ζq(y) p(y) p(y), leading to

(32) Ψ (f ) (ζ) = Y a ζ (y)ν(dy) .
Then, for all y ∈ Y,

|a ζ (y)| sup θ∈T f q(θ, y) p(y) p(y) , (33)
which is integrable w.r.t ν(dy) by (A4)-(iii).

Furthermore, recall that for all y ∈ Y,

[I φ (µ n )q](y) = T µ n (dθ)|b µn (θ)| φ q(θ, y) T µ n (dθ)|b µn (θ)| φ .
By applying twice Lemma 17-(iv) with h(θ) = 1 and h(θ) = q(θ, y), we have that for all y ∈ Y,

lim n→∞ [I φ (µ n )q](y) = [I φ (μ)q](y) . (34)
Now, since f is C 1 by (A2), we obtain from Lemma 17-(i) and (34) respectively that for all y ∈ Y, lim n→∞ a µn (y) = a μ(y) and lim n→∞ a I φ (µn) (y) = a I φ (μ) (y). Combining with [START_REF] Lindsay | Efficiency versus robustness: The case for minimum hellinger distance and related methods[END_REF] and [START_REF] Hellinger | Neue begründung der theorie quadratischer formen von unendlichvielen veränderlichen[END_REF] we can thus apply the Dominated Convergence Theorem to obtain

lim n→∞ Ψ (f ) (µ n ) = Ψ (f ) (μ) (35) and lim n→∞ Ψ (f ) (µ n+1 ) = lim n→∞ Ψ (f ) (I φ (µ n )) = Ψ (f ) (I φ (μ)) . (36)
Finally, [START_REF] Douc | Minimum variance importance sampling via population monte carlo[END_REF] and [START_REF] Royden | Real Analysis (4th Edtion)[END_REF] together yield Ψ (f ) (μ) = Ψ (f ) •I φ (μ), which in turn implies that μ is a fixed point of I φ according to Theorem 1-(ii). Proof of (ii). We prove (ii) by contradiction. Suppose that µ n ⇒ μ, where μ is a fixed point of I φ that satisfies

Ψ (f ) (μ) > inf ζ∈M 1,µ (T) Ψ (f ) (ζ) .
Then, there exists µ ∈ M 1,µ (T) such that Ψ (f ) (μ) > Ψ (f ) (µ ). Now for all n ∈ N, set

B n = θ ∈ T :
µ m (B n 0 ) = Bn 0 µ m-1 (dθ) g µ m-1 (θ) µ m-1 (g µ m-1 ) δµ m-1 (B n 0 ) .
By induction on m we get that, for all m n, µ m (B n 0 ) δ m-n 0 µ n 0 (B n 0 ). This contradicts the previously obtain facts that δ > 1 and µ n 0 (B n 0 ) > 0. Therefore we get a contradiction and the proof is concluded.

A.

A2). Let µ ∈ M 1 (T), φ ∈ R be such that µ(|b µ |) < ∞ and (37) T µ(dθ)E q(θ, Y 1 ) µq(Y 1 ) f µq(Y 1 ) p(Y 1 ) φ < ∞ . Then, (38) lim
K→∞ µ(|b µ,K | φ ) = µ(|b µ | φ ) , P -a.s.
Proof. Set g(θ, y) = q(θ,y) µq(y) |f ( µq(y) p(y) )| and note that E[g(θ, Y 1 )] = |b µ (θ)| since f is of constant sign.

(i) We start with the case φ / ∈ [0, 1]. Our goal is to apply Lemma 22, which is a generalized version of the Dominated Convergence Theorem. To do so, first note that |b µ | φ is positive and combining with the convexity of the mapping u → u φ , we have for all K ∈ N and for all θ ∈ T,

0 |b µ,K (θ)| φ K -1 K k=1 [g(θ, Y k )] φ . (39) Since µ(|b µ |) < ∞, the LLN for µ-almost all θ ∈ T yields lim K→∞ b µ,K (θ) = b µ (θ) . ( 40

  D) := Y p(D, y)ν(dy) is marginal likelihood or model evidence (whose value is unknown). Now denoting by P and µQ the probability measures on (Y, Y) with corresponding associated density p(•|D) and µq with respect to ν and setting p = p(•|D) in Ψ (f ) (µ), the optimisation problem defined by (3) can be rewritten as arginf µ∈M D f (µQ, P) .

Remark 12 .

 12 The Bayesian framework described in Section 2 corresponds to the case p = p(•|D), p = p(D, •) and Z = p(D). Note that in this particular configuration, Lemma 11 allows to upper and lower bound the marginal likelihood [20, Theorem 1].

Algorithm 3 :

 3 α-Approximate f -EI(φ) one-step transition 1. Sampling step : Draw independently Y1, ..., YK ∼ µq 2. Expectation step : bµ,K

Algorithm 4 :

 4 Mixture α-Approximate f -EI(φ) Input: p * : unnormalized version of the density p, Q: Markov transition kernel, K: number of samples, ΘJ = {θ1, ..., θJ } ⊂ T: parameter set. Output: Optimised weights λ. Set λ = [ 1 J , ..., 1 J ]. while the α-bound has not converged do Sampling step : Draw independently K samples Y1, ..., YK from µ λ q.

Figure 1 .

 1 Figure 1. Impact of the hyperparameter φ. Here α = -2 with corresponding range (0, 0.5] for φ and the grey dotted line corresponds to the exact density p.

Algorithm 5 :

 5 Complete Exploitation-Exploration Algorithm Input: p * : unnormalized version of the density p, α: α-divergence parameter, φ: tuning parameter, q0: initial sampler, Q, Qr: Markov transition kernels, K: number of samples, J: dimension of parameter set, (ri)i: rate policy Output: Optimised weights λ and parameter set ΘJ . Draw θ1, ..., θJ from q0. Set i = 0. while the α-bound has not converged do Exploitation step : Set ΘJ = {θ1, ..., θJ }. Perform Mixture α-Approximate f -EI(φ) and obtain λ and ξ (α)

  Toy example revisited. The target p * now corresponds to a mixture of two d-dimensional Gaussian densities multiplied by a positive constant Z such that p * (y) = Z × [0.5N (y; -su d , I d ) + 0.5N (y; su d , I d )]

Figure 2 .

 2 Figure 2. Impact of the dimension d and of the parameter set dimension J. Plotted is the α-bound computed for each pair of (J, d) over 100 replicates.

  (ii) For all ζ ∈ M 1 (T), the function θ → |b ζ (θ)| is continuous. Furthermore for all θ ∈ T, |b µn (θ)| tends to |b μ(θ)| as n → ∞. (iii) There exist 0 < m -< m + < ∞ such that, for all ζ ∈ M 1 (T) and θ ∈ T, |b ζ (θ)| ∈ [m -, m + ].(iv) For all continuous, positive and bounded function h, lim n→∞ T µ n (dθ)|b µn (θ)| φ h(θ) = T μ(dθ)|b μ(θ)| φ h(θ) .

  |b ζ (θ)| = Y a ζ (θ, y)ν(dy) , where we set for all (θ, y) ∈ T×Y, a ζ (θ, y) = q(θ, y) f ζq(y) p(y) (the absolute value can be put inside the integral since f is of constant sign by (A2)). The continuity of |b ζ | follows from the Dominated Convergence Theorem, since for all y ∈ Y, the function θ → a ζ (θ, y) is continuous on T by (A4)-(ii) and for all (θ, y) ∈ T × Y, we have

n→∞a

  µn (θ, y) = a μ(θ, y). The bound (26) and (A4)-(iv) provide a domination criterion and we get that |b µn (θ)| tends to |b μ(θ)| as n → ∞, which concludes the proof of (ii). Proof of (iii). For all (θ, ζ) ∈ T × M 1 (T), we have |b ζ (θ)| ∈ [m -, m + ] where

  Using (ii) the function θ → |b μ(θ)| φ h(θ) is continuous, and, since T is compact,µ n ⇒ μ gives that lim n→∞ T µ n (dθ)|b μ(θ)| φ h(θ) = T μ(dθ)|b μ(θ)| φ h(θ) . (28)Next we show thatlim n→∞ T µ n (dθ) |b µn (θ)| φ -|b μ(θ)| φ h(θ) = 0 (29) Using (iii), since u → u φ is Lipschitz on [m -, m + ], there exists a constant C such that µ n |b µn | φ -|b μ| φ h C sup θ∈T h(θ) T µ n (dθ) ||b µn (θ)| -|b μ(θ)|| = C sup θ∈T h(θ) Y |a n (y)|ν(dy)where a n (y) := µ n q(y) f µnq(y) p(y)-f μq(y) p(y)

  m n {g µm (θ) > δµ m (g µm )} , where for all ζ ∈ M 1 (T), g ζ := |b ζ | φ . There exists, according to Lemma 20, for a well chosen δ > 1, a sufficiently large n 0 such that µ (B n 0 ) > 0. Furthermore µ ≈ µ by definition, where ζ ≈ µ if and only if for all A ∈ T : ζ(A) > 0 is equivalent to µ(A) > 0. Since 0 < |b µ (θ)| φ < ∞ for µ-almost all θ ∈ T and dµ 1 dµ ∝ |b µ | φ , we also have µ 1 ≈ µ. Then by induction, µ n ≈ µ for all n ∈ N. Finally, µ n 0 (B n 0 ) > 0. Moreover, for all θ ∈ B n 0 and all m > n 0 , gµ m (θ) µm(gµ m ) > δ and consequently

)

  Now applying successively (a) the LLN for µ-almost all θ ∈ T (as stated inLemma 23), which is valid under (37), (b) Fubini's Theorem and (c) again the LLN (41) T µ(dθ) lim K→∞ K -1 K k=1 {g(θ, Y k )} φ (a) = T µ(dθ)E {g(θ, Y 1 )} φ •, Y k )] φ < ∞Combining with (39) and (40), we apply Lemma 22 and obtainµ |b µ | φ = µ lim K→∞ |b µ,K | φ = lim K→∞ µ(|b µ,K | φ ) . (ii) We now turn to the case φ ∈ (0, 1]. Let M > 0. Since T µ(dθ) K -1 K k=1 g(θ, Y k )1 {g(θ,Y k ) M } φ µ(|b µ,K | φ ) ,the LLN for µ-almost all θ ∈ T (Lemma 23) and the Dominated Convergence Theorem yields(42) T µ(dθ) E[g(θ, Y 1 )1 {g(θ,Y 1 ) M } ] φ lim inf K→∞ µ(|b µ,K | φ ) .Using now (u + v) φ u φ + v φ and then Jensen's inequality for the concave mapping u → u φ ,µ(|b µ,K | φ ) , Y k )1 {g(θ,Y k ) M } φ , Y k )1 {g(θ,Y k )>M } φBy invoking the LLN for µ-almost all θ ∈ T (Lemma 23) and the Dominated Convergence Theorem for the first term of the rhs and the LLN combined with Fubini for the second term, we getlim sup K→∞ µ(|b µ,K | φ ) T µ(dθ) E[g(θ, Y 1 )1 {g(θ,Y 1 ) M } ] φ + T µ(dθ)E[g(θ, Y 1 )1 {g(θ,Y 1 )>M } ] φLetting M go to infinity both in this inequality and in (42) completes the proof of (38).

Lemma 23 . 1 ZFigure 3 .

 2313 Figure 3. Dimension 1: Impact of the parameter φ. Here α = 0.5 with corresponding range (0, 1] for φ.

Figure 4 .

 4 Figure 4. Dimension 1: Impact of the parameter φ. Here α = 1.2 with corresponding range (-5, 0) for φ.

Table 1

 1 Special cases in the f -divergence family

  3. Lemma 21 : statement and proof. Recall that Y 1 , Y 2 , ... are i.i.d random variables with common density µq w.r.t ν, defined on the same probability space (Ω, F, P) and we denote by E the associated expectation operator.

	Lemma 21. Assume (A1) and (

K. DAUDEL ET AL.

APPENDIX A

A.1. Proof of Theorem 2. In Proposition 5, the difference Ψ (f ) (ζ) -Ψ (f ) (µ) is split into two terms

where g = dζ/dµ. Moreover, Proposition 5 states that A(µ, ζ) is always non-positive.

It turns out that the second term is minimal over all positive probability densities g when it is proportional to |b µ | 1/(1-) , as we show in Lemma 14 below.

Lemma 14. For any positive probability density g w.r.t µ, we have

. Thus Jensen's inequality yields, for any positive probability density g w.r.t. µ,

which finishes the proof of the inequality. The next statement follows from the case of equality in Jensen's inequality: g must be proportional to |b µ | 1/(1-) .

The next lemma shows that this choice leads to a non-positive second term, thus implying that Ψ (f ) (ζ) Ψ (f ) (µ).

Lemma 15. Assume (A1), (A2) and (A3). Then φ = 1/(1 -) satisfies [START_REF] Hoffman | Stochastic variational inference[END_REF]

Proof. We apply [START_REF] Dehaene | Expectation Propagation in the large-data limit[END_REF] with g = 1 and get that

Then [START_REF] Hoffman | Stochastic variational inference[END_REF] can be readily checked with φ = 1/(1 -). Furthermore using µ(|b µ |) < ∞ when φ < 0 and (A1) combined with (A2) for φ > 0, we obtain µ(|b µ | φ ) > 0, which concludes the proof.

and we obtain

which proves [START_REF] Del | On asymptotic properties of informationtheoretic divergences[END_REF].

We now prove [START_REF] Sason | On f-divergences: Integral representations, local behavior, and inequalities[END_REF] in the case where f is non-increasing. First note that for all δ > 1, (δ -1)μ(b μ) 0. Let us define A δ = {b μ < δ μ(b μ)} and show that µ (A δ ) > 0 for some δ > 1. To do so, we proceed by contradiction. Suppose that µ (A δ ) = 0 for all δ > 1, so that

Letting δ ↓ 1, we obtain a contradiction, which finishes the proof. The alternative case where f is non-decreasing is obtained similarly by taking δ ↑ 1.

Lemma 19. Assume (A1) and (A2). Let μ ∈ M 1 (T) be a fixed point of I φ , let g μ := |b μ| φ and φ ∈ R . Let µ ∈ M 1 (T) and assume that there exists µ ∈ M 1,µ (T) such that Ψ (f ) (μ) > Ψ (f ) (µ ). Then, there exists δ > 1 such that µ (g μ > δ μ(g μ)) > 0 in the cases (a) and (b) of Theorem 3.

Proof. Note that (5) holds for any φ and ζ (in particular ζ = μ) by Lemma 17-(iii). As μ is a fixed point of I φ , g μ is μ-almost all constant. Consequently, μ(|b μ| φ ) 1/φ = μ(|b μ|). We separate the two cases f non-increasing and f non-decreasing:

(i) Let f be non-increasing and φ > 0. Then, |b μ| = -b μ and u → u 1/φ is increasing. For all δ > 1, δ := δ 1/φ > 1 and

(ii) Let f be non-decreasing and φ < 0. Then, |b μ| = b μ and u → u 1/φ is decreasing. For all δ > 1, δ := δ 1/φ < 1 and

We conclude by applying Lemma 18 in the two separated cases.

Lemma 20. Assume (A1), (A2) and (A4). Let φ ∈ R and denote

and define the sequence (µ n ) n∈N according to (4). Suppose that µ n ⇒ μ for some fixed point μ ∈ M 1 (T) of I φ . Further assume there exists µ ∈ M 1,µ (T) such that

Then, there exist δ > 1 and n ∈ N * such that

in the cases (a) and (b) of Theorem 3.

Proof. First note that the sequence (µ n ) n∈N is well-defined for any φ ∈ R by Lemma 17-(iii), which implies µ n (g µn ) > 0 for all n. We further have that

Furthermore, applying (ii) and (iv) in Lemma 17, we have, for all θ ∈ T, lim n→∞ g µn (θ) = g μ(θ) and lim n→∞ µ n (g µm ) = μ(g μ). Hence, for all θ ∈ T,

The proof is concluded by applying Lemma 19.

Proof of Theorem 3. Assume (A1), (A2) and (A4). Lemma 17-(iii) is exactly the first result we want to obtain, that is: for all

APPENDIX B: TECHNICAL RESULTS

B.1. General Dominated Convergence Theorem. We state and prove a generalized version of the Dominated Convergence Theorem, adapted from [START_REF] Royden | Real Analysis (4th Edtion)[END_REF]Theorem 19]. We provide here a full proof for the sake of completeness.