# Monotonic Alpha-divergence Minimisation for Variational Inference

Kamélia Daudel, Randal Douc, François Roueff

<figure><img src="image_1.png" /></figure>

## To cite this version:

Kamélia Daudel, Randal Douc, François Roueff. Monotonic Alpha-divergence Minimisation for Variational Inference. Journal of Machine Learning Research, 2023, 24 (62), pp.1-76. hal-03164338v2

HAL Id: hal-03164338

https://telecom-paris.hal.science/hal-03164338v2

Submitted on 27 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Monotonic Alpha-divergence Minimisation for Variational Inference

**Kamélia Daudel**¹,²

1: LTCI, Télécom Paris

Institut Polytechnique de Paris, France

2: Department of Statistics, University of Oxford

Oxford OX1 3TG, United Kingdom

KAMELIA.DAUDEL@GMAIL.COM

**Randal Douc**

SAMOVAR, Télécom SudParis

Institut Polytechnique de Paris, France

RANDAL.DOUC@TELECOM-SUDPARIS.EU

**François Roueff**

LTCI, Télécom Paris

Institut Polytechnique de Paris, France

FRANCOIS.ROUEFF@TELECOM-PARIS.FR

**Editor:** Shakir Mohamed

## Abstract

In this paper, we introduce a novel family of iterative algorithms which carry out α-divergence minimisation in a Variational Inference context. They do so by ensuring a systematic decrease at each step in the α-divergence between the variational and the posterior distributions. In its most general form, the variational distribution is a mixture model and our framework allows us to simultaneously optimise the weights and components parameters of this mixture model. Our approach permits us to build on various methods previously proposed for α-divergence minimisation such as Gradient or Power Descent schemes and we also shed a new light on an integrated Expectation Maximization algorithm. Lastly, we provide empirical evidence that our methodology yields improved results on several multimodal target distributions and on a real data example.

**Keywords:** Variational Inference, Kullback-Leibler, Alpha-Divergence, Mixture Models, Bayesian Inference

## 1. Introduction

Bayesian inference tasks often induce intractable and hard-to-compute posterior densities which need to be approximated. Among the class of approximating methods, Variational inference methods (for example Variational Bayes Jordan et al., 1999; James Beal, 2003) have attracted a lot of attention as they have empirically been shown to be widely applicable to many high-dimensional machine-learning problems (Hoffman et al., 2013; Ranganath et al., 2014; Kingma and Welling, 2014).

These optimisation-based methods introduce a simpler variational family $\mathcal{Q}$ and find the best approximation to the unknown posterior density among this family in terms of a
certain measure of dissimilarity, the most common choice of measure of dissimilarity being
the exclusive Kullback-Leibler divergence (Blei et al., 2017; Zhang et al., 2019).

However, the exclusive Kullback-Leibler divergence is known to have some drawbacks.
Indeed, its zero-forcing behavior is responsible for returning variational approximations
with light tails that severely underestimate the posterior covariance and that are unable to
capture multimodality within the posterior density (Minka, 2005; Jerfel et al., 2021; Prangle,
2021). This is especially inconvenient when the variational family $\mathcal{Q}$ does not exactly match
the posterior density (Yao et al., 2018; Campbell and Li, 2019) and even more so if one
wants to appeal to Importance Sampling methods to approximate integrals of interest in a
Bayesian Inference setting (Jerfel et al., 2021; Prangle, 2021).

To avoid this hurdle, advances in Variational Inference turned to more general families of
divergences such as the $\alpha$-divergence (Zhu and Rohwer, 1995b,a) and Rényi's $\alpha$-divergence
(Rényi, 1961; van Erven and Harremoes, 2014). These families of divergences both recover
the exclusive Kullback-Leibler when $\alpha \to 1$ and thanks to the hyperparameter $\alpha$, they
provide a more flexible framework that can bypass the difficulties associated to the exclusive
Kullback-Leibler divergence by choosing $\alpha < 1$ (Minka, 2005). For this reason, they have
notably been used in Minka (2004); Hernandez-Lobato et al. (2016); Li and Turner (2016);
Dieng et al. (2017); Wang et al. (2018); Daudel et al. (2021); Daudel and Douc (2021).

In the spirit of Variational Inference methods based on the $\alpha$-divergence, we propose in
this paper to build a framework for $\alpha$-divergence minimisation in a Variational Inference
context. The particularity of our work is that our algorithms will ensure a monotonic decrease
at each step in the $\alpha$-divergence between the variational and the posterior distributions. In
addition, our work will apply to variational families as wide as the class of mixture models.
The paper is then organised as follows:

• In Section 2, we introduce some notation and state the optimisation problem we aim at solving in terms of the posterior density, the variational density $q \in \mathcal{Q}$ and the $\alpha$-divergence.

• In Section 3, we consider the typical Variational Inference case where $q$ belongs to a parametric family. In this particular case, we state in Theorem 1 conditions which ensure a systematic decrease in the $\alpha$-divergence at each step for all $\alpha \in [0, 1)$. We then show in Corollary 1 that these conditions are satisfied for a well-chosen iterative scheme and we call the resulting approach the *maximisation approach*. This approach is particularly convenient, a fact that we illustrate over Example 1 and Lemma 1. Furthermore, we derive in Corollary 2 additional iterative schemes satisfying the conditions of Theorem 1 under the name *gradient-based approach*, which we then use to underline the links between our approach and Gradient Descent schemes for $\alpha$-divergence and Rényi's $\alpha$-divergence minimisation.

• In Section 4, we extend the results from Section 3 to the more general case of mixture
models. We derive in Theorems 2 and 3 conditions to simultaneously optimise the weights
and the components parameters of a given mixture model, all the while maintaining the
systematic decrease in the $\alpha$-divergence initially enjoyed in Section 3. These conditions are
then met in Corollary 4 and 5, which respectively generalise the maximisation approach of
Corollary 1 and the gradient-based approach of Corollary 2 to the case of mixture models.

• Section 5 is devoted to related work. We explain in more detail the links between our
framework and existing Variational Inference methods for $\alpha$-divergence minimisation. We
also connect our approach to the Power Descent algorithm from Daudel et al. (2021) and provide in Proposition 2 additional monotonicity results which go beyond the case $\alpha \in [0, 1)$. In addition, we obtain that an integrated Expectation Maximization algorithm introduced in Cappé et al. (2008) can be recovered as a special case of our framework.

* In Section 6, we state generic results that solve our maximisation and gradient-based approaches when the variational family is based on the exponential family. More specifically, we introduce in Section 6.1 additional notation for the exponential family and we recall several of its useful properties. Section 6.2 and Section 6.3 then focus on the maximisation approach and gradient-based approach respectively and notably provide the theoretical justification behind some special cases mentioned earlier in the paper, such as Gaussian densities in Example 1. Lastly, Section 6.4 extends once again the maximisation approach in order to incorporate variational families that were not previously included in our framework, yet draw on the exponential family (such as multivariate Student's t densities in Example 5).

* In Section 7, we focus on the important case of Gaussian Mixture Models (GMMs) and we discuss practical implementations of our algorithms for GMMs optimisation.

* Finally, we show in Section 8 that our enhanced framework has empirical benefits that permit us to improve on the existing algorithms presented in Section 5 when we consider a variety of multimodal targets as well as a real data example.

## 2. Notation and Optimisation Problem

Let $(\mathcal{Y}, \mathcal{Y}, \nu)$ be a measured space, where $\nu$ is a $\sigma$-finite measure on $(\mathcal{Y}, \mathcal{Y})$. Assume that we have access to some observed variables $\mathcal{D}$ generated from a probabilistic model $p(\cdot|y)$ parameterised by a hidden random variable $y \in \mathcal{Y}$ that is drawn from a certain prior $p_0$. The posterior density of the latent variable $y$ given the data $\mathcal{D}$ is then given by

$$p(y|\mathcal{D}) = \frac{p(y, \mathcal{D})}{p(\mathcal{D})} = \frac{p_0(y)p(\mathcal{D}|y)}{p(\mathcal{D})},$$

where $p(\mathcal{D}) = \int_{\mathcal{Y}} p_0(y)p(\mathcal{D}|y)\nu(dy)$ is called the marginal likelihood or model evidence. For many interesting choices of models, sampling directly from the posterior density is impossible and the marginal likelihood is also unknown or too costly to compute.

To address this difficulty, Variational Inference methods approximate the posterior density by a simpler probability density belonging to a given variational family $\mathcal{Q}$ (from which it is easy to sample from). They select the best probability density in $\mathcal{Q}$ by solving an optimisation problem that involves a certain measure of dissimilarity, which is chosen to be the $\alpha$-divergence in this paper due to its advantages compared to the more traditional exclusive Kullback-Leibler divergence (see for example Minka, 2005).

More precisely, let us denote by $\mathbb{P}$ the probability measure on $(\mathcal{Y}, \mathcal{Y})$ with corresponding density $p(\cdot|\mathcal{D})$ with respect to $\nu$. As for the variational family, let us denote by $\mathbb{Q}$ the probability measure on $(\mathcal{Y}, \mathcal{Y})$ with associated density $q \in \mathbb{Q}$ with respect to $\nu$. We let $f_\alpha$ be the convex function on $(0, +\infty)$ defined by $f_0(u) = -\log(u)$, $f_1(u) = u \log u$ and $f_\alpha(u) = \frac{1}{\alpha(\alpha-1)}[u^\alpha - 1]$ for all $\alpha \in \mathbb{R} \setminus \{0, 1\}$. Then, assuming that both $\mathbb{Q}$ and $\mathbb{P}$ are
absolutely continuous w.r.t. $\nu$, the $\alpha$-divergence between $\mathbb{Q}$ and $\mathbb{P}$ (extended by continuity to the cases $\alpha = 0$ and $\alpha = 1$ as for example done in Cichocki and Amari, 2010) is given by

$$D_{\alpha}(\mathbb{Q}||\mathbb{P}) = \int_{\mathcal{Y}} f_{\alpha}\left(\frac{q(y)}{p(y|\mathcal{D})}\right) p(y|\mathcal{D})\nu(dy), \quad (1)$$

and the Variational Inference optimisation problem we aim at solving is $\inf_{q \in \mathcal{Q}} D_{\alpha}(\mathbb{Q}||\mathbb{P})$. Notably, it can easily be proven that this optimisation problem is equivalent to solving

$$\inf_{q \in \mathcal{Q}} \Psi_{\alpha}(q; p) \quad \text{with} \quad p(y) = p(y, \mathcal{D}) \quad \text{for all } y \in \mathcal{Y}, \qquad (2)$$

where, for all measurable non-negative function $p$ on $(\mathcal{Y}, \mathcal{Y})$ and for all $q \in \mathcal{Q}$, we have set

$$\Psi_{\alpha}(q; p) = \int_{\mathcal{Y}} f_{\alpha}\left(\frac{q(y)}{p(y)}\right) p(y)\nu(dy).$$

As the unknown constant $p(\mathcal{D})$ does not appear anymore in the optimisation problem (2), this formulation is often preferred. Therefore, we consider the general optimisation problem

$$\inf_{q \in \mathcal{Q}} \Psi_{\alpha}(q; p), \qquad (3)$$

where $p$ is any measurable positive function on $(\mathcal{Y}, \mathcal{Y})$. Note that we may drop the dependency on $p$ in $\Psi_{\alpha}$ for notational ease and when no ambiguity occurs.

At this stage, we are left with the choice of the variational family $\mathcal{Q}$ appearing in the optimisation problem (3). The natural idea in Variational Inference and the starting point of our approach is then to work within a parametric family : letting $(\mathcal{T}, \mathcal{T})$ be a measurable space, $K : (\theta, A) \mapsto \int_{A} k(\theta, y)\nu(dy)$ be a Markov transition kernel on $\mathcal{T} \times \mathcal{Y}$ with kernel density $k$ defined on $\mathcal{T} \times \mathcal{Y}$, we consider a parametric family of the form

$$\mathcal{Q} = \{q : y \mapsto k(\theta, y) : \theta \in \mathcal{T}\}.$$

### 3. An Iterative Algorithm for Optimising $\Psi_{\alpha}(k(\theta, \cdot); p)$

In this section, our goal is to define iterative procedures which optimise $\Psi_{\alpha}(k(\theta, \cdot); p)$ with respect to $\theta$ and which are such that they ensure a systematic decrease in $\Psi_{\alpha}$ at each step. For this purpose, we start by introducing some mild conditions on $k$, $p$ and $\nu$ that will be used throughout the paper.

(A1) The kernel density $k$ on $\mathcal{T} \times \mathcal{Y}$, the function $p$ on $\mathcal{Y}$ and the $\sigma$-finite measure $\nu$ on $(\mathcal{Y}, \mathcal{Y})$ satisfy, for all $(\theta, y) \in \mathcal{T} \times \mathcal{Y}$, $k(\theta, y) > 0$, $p(y) \ge 0$ and $0 < \int_{\mathcal{Y}} p(y)\nu(dy) < \infty$.

From this point onwards, and as announced in the previous section, we will drop the dependency on $p$ in $\Psi_{\alpha}$ for notational convenience. Let us now construct a sequence $(\theta_n)_{n \ge 1}$ valued in $\mathcal{T}$ such that the sequence $(\Psi_{\alpha}(k(\theta_n, \cdot)))_{n \ge 1}$ is decreasing. The core idea of our approach will rely on the following proposition.
**Proposition 1** Assume (A1). For all $\alpha \in [0, 1)$ and all $\theta, \theta' \in \mathbb{T}$ such that $\Psi_{\alpha}(k(\theta', \cdot)) < \infty$, it holds that

$$ \Psi_{\alpha}(k(\theta, \cdot)) \le \int_{\mathbb{Y}} \frac{k(\theta', y)^{\alpha} p(y)^{1-\alpha}}{\alpha - 1} \log \left( \frac{k(\theta, y)}{k(\theta', y)} \right) \nu(dy) + \Psi_{\alpha}(k(\theta', \cdot)). \quad (4) $$

Moreover, if $\alpha = 0$, this inequality is an equality and, if $\alpha \in (0, 1)$, equality in (4) is equivalent to having $k(\theta, y) = k(\theta', y)$ for $\nu$-a.e. $y \in \{p > 0\}$.

**Proof** We treat the two cases $\alpha = 0$ and $\alpha \in (0, 1)$ separately.

(a) Case $\alpha = 0$ with $f_0(u) = -\log(u)$ for all $u > 0$. This case is immediate since

$$ \Psi_0(k(\theta, \cdot)) = - \int_{\mathbb{Y}} p(y) \log \left( \frac{k(\theta, y)}{k(\theta', y)} \right) \nu(dy) + \Psi_0(k(\theta', \cdot)). $$

(b) Case $\alpha \in (0, 1)$ with $f_{\alpha}(u) = \frac{1}{\alpha(\alpha-1)}[u^{\alpha} - 1]$ for all $u > 0$. We have that

$$ \begin{aligned} \Psi_{\alpha}(k(\theta, \cdot)) &= \int_{\mathbb{Y}} \left[ \left( \frac{k(\theta, y)}{p(y)} \right)^{\alpha} - 1 \right] p(y) \nu(dy) \\ &= \int_{\mathbb{Y}} \left( \frac{k(\theta', y)}{p(y)} \right)^{\alpha} \frac{\left[ \left( \frac{k(\theta, y)}{k(\theta', y)} \right)^{\alpha} - 1 \right]}{\alpha(\alpha-1)} p(y) \nu(dy) + \Psi_{\alpha}(k(\theta', \cdot)) \end{aligned} $$

Furthermore, using that $\log(u^{\alpha}) \le u^{\alpha} - 1$ for all $u > 0$ with equality only if $u = 1$, and since $\alpha \in (0, 1)$, this inequality becomes

$$ \frac{u^{\alpha} - 1}{\alpha(\alpha - 1)} \le \frac{\log(u^{\alpha})}{\alpha(\alpha - 1)} = \frac{\log(u)}{\alpha - 1}. $$

Applying this with $u = k(\theta, y)/k(\theta', y)$ in the previous expression, we get that

$$ \Psi_{\alpha}(k(\theta, \cdot)) \le \int_{\mathbb{Y}} \frac{k(\theta', y)^{\alpha} p(y)^{1-\alpha}}{\alpha - 1} \log \left( \frac{k(\theta, y)}{k(\theta', y)} \right) \nu(dy) + \Psi_{\alpha}(k(\theta', \cdot)) $$

which is exactly (4). Moreover, the equality holds if $k(\theta, y)/k(\theta', y) = 1$ for $\nu$-a.e. $y$ such that $k(\theta', y)^{\alpha} p(y)^{1-\alpha} > 0$, which, by (A1), only happens if $p(y) > 0$.

<figure><img src="image_3.png" /></figure>

This result then allows us to deduce Theorem 1 below.

**Theorem 1** Assume (A1). Let $\alpha \in [0, 1)$ and starting from an initial $\theta_1 \in \mathbb{T}$, let $(\theta_n)_{n \ge 1}$ be defined iteratively such that for all $n \ge 1$,

$$ \int_{\mathbb{Y}} \frac{k(\theta_n, y)^{\alpha} p(y)^{1-\alpha}}{\alpha - 1} \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_n, y)} \right) \nu(dy) \le 0. \quad (5) $$

Further assume that $\Psi_{\alpha}(k(\theta_1, \cdot)) < \infty$. Then, $\Psi_{\alpha}(k(\theta_{n+1}, \cdot)) \le \Psi_{\alpha}(k(\theta_n, \cdot))$ for all $n \ge 1$ with an equality occurring only if (5) is also an equality.
**Proof** Applying Proposition 1 iteratively for $n = 1, 2, \dots$ with $\theta = \theta_{n+1}$ and $\theta' = \theta_n$ in (4) combined with (5), we get that $\Psi_\alpha(k(\theta_n, \cdot)) < \infty$ and

$$ \Psi_{\alpha}(k(\theta_{n+1}, \cdot)) \le \int_{Y} \frac{k(\theta_n, y)^{\alpha} p(y)^{1-\alpha}}{\alpha - 1} \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_n, y)} \right) \nu(dy) + \Psi_{\alpha}(k(\theta_n, \cdot)) \le \Psi_{\alpha}(k(\theta_n, \cdot)). $$

Suppose now that $\Psi_\alpha(k(\theta_{n+1}, \cdot)) = \Psi_\alpha(k(\theta_n, \cdot))$ for some $n \ge 1$. Then we have equalities in the previous display and thus an equality in (5) and in (4) with $\theta = \theta_{n+1}$ and $\theta' = \theta_n$. ■

At this point, we seek to find iterative schemes satisfying (5). To do so, for all $\alpha \in [0, 1)$, all $n \ge 1$ and all $y \in Y$, we introduce the notation:

$$ \varphi_n^{(\alpha)}(y) = k(\theta_n, y)^{\alpha} p(y)^{1-\alpha} \quad (6) $$

$$ \check{\varphi}_n^{(\alpha)} = \frac{\varphi_n^{(\alpha)}}{\int \varphi_n^{(\alpha)} d\nu}. $$

Note that under (A1), the normalisation in $\check{\varphi}_n^{(\alpha)}$ satisfies $0 < \int \varphi_n^{(\alpha)} d\nu < \infty$, where the right inequality in particular follows from Jensen's inequality applied to the strictly concave function $u \mapsto u^{1-\alpha}$. We then state our first corollary.

**Corollary 1 (Maximisation approach)** Assume (A1). Let $\alpha \in [0, 1)$ and let $(b_n)_{n \ge 1}$ be a non-negative sequence. Starting from an initial $\theta_1 \in T$, let $(\theta_n)_{n \ge 1}$ be defined iteratively as follows

$$ \theta_{n+1} = \underset{\theta \in T}{\operatorname{argmax}} \int_{Y} \left[ \varphi_n^{(\alpha)}(y) + b_n k(\theta_n, y) \right] \log \left( \frac{k(\theta, y)}{k(\theta_n, y)} \right) \nu(dy), \quad n \ge 1, \quad (7) $$

where we assume that this argmax is uniquely defined at each step. Then (5) holds at all times $n \ge 1$ and we can apply Theorem 1.

**Proof** Under (A1) and using that $\alpha \in [0, 1)$, we can rewrite (5) as

$$ \int_{Y} \varphi_{n}^{(\alpha)}(y) \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_{n}, y)} \right) \nu(dy) \ge 0. $$

Since $\int_Y k(\theta_n, y) \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_n, y)} \right) \nu(dy) = -D_1(K(\theta_n, \cdot) || K(\theta_{n+1}, \cdot)) \le 0$, we also obtain that

$$ \int_{Y} \varphi_{n}^{(\alpha)}(y) \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_{n}, y)} \right) \nu(dy) \ge \int_{Y} \left[ \varphi_{n}^{(\alpha)}(y) + b_{n} k(\theta_{n}, y) \right] \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_{n}, y)} \right) \nu(dy). $$

Thus, (5) holds at all times $n$ if the following condition is satisfied:

$$ \int_{Y} \left[ \varphi_{n}^{(\alpha)}(y) + b_{n} k(\theta_{n}, y) \right] \log \left( \frac{k(\theta_{n+1}, y)}{k(\theta_{n}, y)} \right) \nu(dy) \ge 0, \quad n \ge 1, $$

which is implied by the definition of $\theta_{n+1}$ in (7).

We now make three comments regarding Corollary 1.
• The r.h.s. of (7) can be simplified, yielding the equivalent iterative scheme

$$ \theta_{n+1} = \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_{\mathcal{Y}} \left[ \varphi_n^{(\alpha)}(y) + b_n k(\theta_n, y) \right] \log k(\theta, y) \nu(dy), \quad n \ge 1, $$

provided that $\int_{\mathcal{Y}} [\varphi_n^{(\alpha)}(y) + b_n k(\theta_n, y)] |\log k(\theta_n, y)| \nu(dy) < \infty$ for all $n \ge 1$.

• While it suffices to find any $\theta_{n+1}$ so that (5) holds to obtain a systematic decrease in $\Psi_\alpha$, defining $\theta_{n+1}$ as in (7) enables us to solve this argmax problem for notable choices of kernel density $k$. A remarkable aspect is indeed that (7) is written as a maximisation problem involving the logarithm of the ratio $k(\theta, y)/k(\theta_n, y)$ whereas $\Psi_\alpha$ is not directly expressed in terms of logarithm of this ratio for $\alpha \notin \{0, 1\}$. As a result, we can use the maximisation approach of Corollary 1 to derive simple update rules for $(\theta_n)_{n \ge 1}$. An example and a lemma are provided thereafter to illustrate this fact and more general results regarding exponential families will follow later in Section 6.

• The sequence $(b_n)_{n \ge 1}$ appearing in (7) is arbitrary and can as a result be chosen by the practitioner. It is responsible for introducing a regularisation term in the argmax problem (7) so that the larger $b_n$ is, the closer $\theta_{n+1}$ is to $\theta_n$.

We next provide a motivating example where the maximisation approach is applicable.

**Example 1 (Maximisation approach for a Gaussian density)** We consider a d-dimensional Gaussian density, in which case $\mathcal{Y} = \mathbb{R}^d$, $\nu$ is the d-dimensional Lebesgue measure and $k(\theta, y) = \mathcal{N}(y; m, \Sigma)$, where $\theta = (m, \Sigma) \in \mathcal{T}$ denotes the mean and covariance matrix of the Gaussian density. We let the parameter set $\mathcal{T}$ include all possible means $m$ in $\mathbb{R}^d$ and all possible positive-definite covariance matrices $\Sigma$. Finally, denoting by $\|\cdot\|$ the Euclidean norm, we assume that the non-negative function $p$ defined on $\mathcal{Y}$ satisfies

$$ 0 < \int_{\mathcal{Y}} \left(1 + \|y\|^{2/(1-\alpha)}\right) p(y) \, dy < \infty. \qquad (8) $$

Then (A1) holds. Furthermore, starting from any $\theta_1 = (m_1, \Sigma_1) \in \mathcal{T}$ so that $\Psi_\alpha(k(\theta_1, \cdot)) < \infty$ and denoting $\theta_n = (m_n, \Sigma_n)$ for all $n \ge 1$, it holds that: for all $n \ge 1$ and all $\gamma_n \in (0, 1]$, there exists $b_n \ge 0$ such that the argmax problem (7) has a unique solution defined by

$$ m_{n+1} = \gamma_n \hat{m}_n + (1 - \gamma_n) m_n, \\ \Sigma_{n+1} = \gamma_n \hat{\Sigma}_n + (1 - \gamma_n) \Sigma_n + \gamma_n (1 - \gamma_n) (\hat{m}_n - m_n) (\hat{m}_n - m_n)^T, \qquad (9) $$

where, using the definition of $\tilde{\varphi}_n^{(\alpha)}$ in (6), we set

$$ \hat{m}_n = \int_{\mathcal{Y}} y \tilde{\varphi}_n^{(\alpha)}(y) \nu(dy), \\ \hat{\Sigma}_n = \int_{\mathcal{Y}} y y^T \tilde{\varphi}_n^{(\alpha)}(y) \nu(dy) - \hat{m}_n \hat{m}_n^T. $$

Here, the detailed derivation of (9) is deferred to Section 6.2 as we focus on interpreting this result in light of Corollary 1 (in particular we will see in Corollary 6 that the case
$\gamma_n = 1$ corresponding to $b_n = 0$ will require an additional non-degenerate condition on $p$). By Corollary 1, the systematic decrease in $(\Psi_\alpha(k(\theta_n, \cdot)))_{n \ge 1}$ holds for any choice of sequence $(\gamma_n)_{n \ge 1}$ valued in $(0, 1]$ in the updates (9). In addition, $\gamma_n$ permits us to build a tradeoff at time $n \ge 1$ between selecting an update close to the current parameter (that is, taking $\gamma_n$ close to zero) and choosing the Gaussian density with exactly the same mean and covariance matrix as $\check{\varphi}_n^{(\alpha)}$ (that is, setting $\gamma_n = 1$). This is a key idea that is linked to the regularisation term appearing in (7) and that we will revisit several times throughout the paper.

The maximisation approach of Corollary 1 is also applicable to the commonly-used mean-field variational family, which approximates the true density by a density with independent components parameterised separately. The following lemma indeed states that the global argmax problem can be separated into component-wise argmax problems in that case.

**Lemma 1 (Maximisation approach for the mean-field family)** A generic member of the mean-field variational family is given by $k(\theta, y) = \prod_{\ell=1}^L k^{(\ell)}(\theta^{(\ell)}, y^{(\ell)})$ with $\theta = (\theta^{(1)}, \dots, \theta^{(L)}) \in \mathbb{T}$ and $y = (y^{(1)}, \dots, y^{(L)}) \in \mathbb{Y}$. Then, starting from $\theta_1 \in \mathbb{T}$ and denoting $\theta_n = (\theta_n^{(1)}, \dots, \theta_n^{(L)})$ for all $n \ge 1$, solving (7) yields the following update formulas:

$$ \theta_{n+1}^{(\ell)} = \underset{\theta^{(\ell)}}{\operatorname{argmax}} \int_{\mathbb{Y}} \left[ \varphi_n^{(\alpha)}(y) + b_n k(\theta_n, y) \right] \log \left( \frac{k^{(\ell)}(\theta^{(\ell)}, y^{(\ell)})}{k^{(\ell)}(\theta_n^{(\ell)}, y^{(\ell)})} \right) \nu(dy), \quad \ell = 1 \dots L, \quad n \ge 1. $$

The maximisation approach is not the only way to satisfy (5). Indeed, this can also be achieved by taking a gradient-based approach and relying on additional smoothness conditions (see Appendix A.1 for the definition of $\beta$-smoothness), as written in Corollary 2.

**Corollary 2 (Gradient-based approach)** Assume (A1). Let $\mathbb{T} \subseteq \mathbb{R}^d$ be a convex set, let $\alpha \in [0, 1)$ and let $(\gamma_n)_{n \ge 1}$ be valued in $(0, 1]$. Starting from an initial $\theta_1 \in \mathbb{T}$, let $(\theta_n)_{n \ge 1}$ be defined iteratively as follows

$$ \theta_{n+1} = \theta_n - \frac{\gamma_n}{\beta_n} \nabla g_n(\theta_n), \quad n \ge 1, \qquad (10) $$

where $(g_n)_{n \ge 1}$ is the sequence of functions defined by: for all $n \ge 1$ and all $\theta \in \mathbb{T}$,

$$ g_n(\theta) = \int_{\mathbb{Y}} \frac{\varphi_n^{(\alpha)}(y)}{\alpha - 1} \log \left( \frac{k(\theta, y)}{k(\theta_n, y)} \right) \nu(dy) \qquad (11) $$

and $g_n$ is assumed to be $\beta_n$-smooth. Then (5) holds for all $n \ge 1$ and we can apply Theorem 1.

**Proof** Since $\gamma_n \in (0, 1]$ and $g_n$ is a $\beta_n$-smooth function by assumption, we can apply Lemma 2 and we obtain that for all $n \ge 1$,

$$ g_n(\theta_n) - g_n\left(\theta_n - \frac{\gamma_n}{\beta_n} \nabla g_n(\theta_n)\right) \ge \frac{\gamma_n}{2\beta_n} \|\nabla g_n(\theta_n)\|^2. $$

Thus, by definition of $\theta_{n+1}$ in (10), we have $0 = g_n(\theta_n) \ge g_n(\theta_{n+1})$, which in turn implies (5) and the proof is concluded. ■
Let us now reflect on the implications of Corollary 2. Under common differentiability assumptions, we can write: for all $n \ge 1$ and all $\theta' \in \mathbb{T}$

$$ \nabla g_n(\theta') = \int_{\mathbb{Y}} \frac{\varphi_n^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta', y)} \nu(dy). $$

Then, (10) becomes

$$ \theta_{n+1} = \theta_n - \frac{\gamma_n}{\beta_n} \int_{\mathbb{Y}} \frac{\varphi_n^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy), \quad n \ge 1. \quad (12) $$

Under this form, the iterative scheme (12) bears similarities with Gradient Descent iterations for $\alpha$-divergence and Rényi's $\alpha$-divergence minimisation. Indeed, given a learning rate policy ($r_n$)$_{n \ge 1}$ and setting $p = p(\cdot, \mathcal{D})$, such Gradient Descent iterations are given respectively by

$$ \theta_{n+1} = \theta_n - r_n \int_{\mathbb{Y}} \frac{\varphi_n^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy), \quad n \ge 1 $$

$$ \theta_{n+1} = \theta_n - r_n \int_{\mathbb{Y}} \frac{\check{\varphi}_n^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy), \quad n \ge 1 $$

(we refer to Appendix A.2 for details regarding how these updates are obtained). Building on this comment, let us now give an example where the conditions on ($g_n$)$_{n \ge 1}$ from Corollary 2 are satisfied and show how Gradient Descent steps (in that case for Rényi's $\alpha$-divergence minimisation) can originate from our gradient-based approach.

**Example 2 (Gradient-based approach for a Gaussian density)** We consider the case of a $d$-dimensional Gaussian density with $k(\theta, y) = N(y; m, \Sigma)$, where this time $\theta = m$. Further assume that $\mathbb{T} \subseteq \mathbb{R}^d$ is a convex subset and that $\Sigma = \sigma^2 \mathbf{I}_d$, where $\sigma^2 > 0$ and $\mathbf{I}_d$ denotes the $d$-dimensional identity matrix. Finally, denoting by $\|\cdot\|$ the Euclidean norm, we assume that the non-negative function $p$ defined on $\mathbb{Y}$ satisfies

$$ 0 < \int_{\mathbb{Y}} \left(1 + \|y\|^{1/(1-\alpha)}\right) p(y) dy < \infty. \quad (13) $$

Then, (A1) holds. Furthermore, denoting $\theta_n = m_n$ and setting $\beta_n = \sigma^{-2}(1-\alpha)^{-1} \int \varphi_n^{(\alpha)} d\nu$ for all $n \ge 1$, the conditions on ($g_n$)$_{n \ge 1}$ from Corollary 2 are satisfied so that the mean of the Gaussian density $N(y; m, \Sigma)$ can be optimised as follows: for all $n \ge 1$ and all $\gamma_n \in [0, 1)$

$$ m_{n+1} = \gamma_n \hat{m}_n + (1 - \gamma_n) m_n, \quad (14) $$

where $\hat{m}_n = \int_{\mathbb{Y}} y \check{\varphi}_n^{(\alpha)}(y) \nu(dy)$ (we refer to Corollary 8 for the detailed derivation). Setting $p = p(\cdot, \mathcal{D})$, the update (14) can notably be seen as a Gradient Descent step for Rényi's $\alpha$-divergence minimisation with a learning rate $r_n = \sigma^2(1-\alpha)\gamma_n$, where $\gamma_n \in [0, 1)$.

We now make two important comments.
* Corollary 2 sheds light on the links between our approach and the more traditional Gradient Descent methodology for optimising objective functions based on the α-divergence in Variational Inference (Li and Turner, 2016; Dieng et al., 2017).

Unlike the usual Gradient Descent methodology, Corollary 2 requires a smoothness condition on $g_n$. Smoothness conditions are tools that are often used to obtain stronger convergence guarantees for Variational Inference algorithms (see for example Buchholz et al., 2018; Alquier and Ridgway, 2020). Yet, it can be difficult to satisfy these smoothness assumptions in practice, even if some results have been derived for specific variational families when using the exclusive Kullback-Leibler divergence as the objective function (Domke, 2020).

To the best of our knowledge, no such results have been proven for α-divergence minimisation. In our case, we obtain that a smoothness condition on $g_n$ translates into a systematic decrease in $\Psi_α$. As we shall see in the forthcoming section, being able to establish monotonicity results on $(\theta_n)_{n \ge 1}$ will come in handy as we try to go beyond the framework considered in Section 3.

* To put things into perspective with the maximisation approach from earlier, observe that the updates on the means in Examples 1 and 2 coincide. The difference between the two examples is due to the fact that the former provides an update for the covariance matrix as well, without sacrificing the monotonicity of the overall algorithm. It is in fact hard to derive a covariance matrix update in the Gaussian case using the gradient-based approach (see Section 6.3 for details). In that sense, the maximisation approach provides an interesting alternative that can bypass this difficult smoothness assumption on $g_n$. We will further delve into this aspect as we reach Section 6.

So far, the log function appearing in Theorem 1 considerably eased the derivation of iterative update formulas for well-chosen kernel densities $k$ that can exploit this log structure (Examples 1-2 and Lemma 1). Yet, these choices of kernel densities can be too restrictive to fully capture complex and multimodal posterior densities. Since working with a larger variational family might lead to more accurate approximations of the posterior density, an idea is then to investigate whether the iterative update formulas from Section 3 can be generalised to the mixture model variational family (for example whether we can extend the updates for a Gaussian kernel from Example 1 to Gaussian Mixture Models).

As we shall see in the next section, further theoretical developments will be required in order to derive valid iterative schemes that optimise both the mixture weights and the mixture components parameters of a given mixture model.

## 4. Extension to Mixture Models

Let us first formally define the class of mixture models we are going to be working with. Given $J \in \mathbb{N}^*$, we introduce the simplex of $\mathbb{R}^J$:

$$ S_J = \left\{ \lambda = (\lambda_1, \dots, \lambda_J) \in \mathbb{R}^J : \forall j \in \{1, \dots, J\}, \lambda_j \ge 0 \text{ and } \sum_{j=1}^J \lambda_j = 1 \right\}, $$
we define $S_J^+ = \{\lambda \in S_J : \forall j \in \{1, \dots, J\}, \lambda_j > 0\}$ and we denote $\Theta = (\theta_1, \dots, \theta_J) \in \mathbb{T}^J$.
We consider the mixture model variational family given by

$$Q = \left\{ q : y \mapsto \mu_{\lambda, \Theta} k(y) := \sum_{j=1}^{J} \lambda_j k(\theta_j, y) : \lambda \in S_J, \Theta \in \mathbb{T}^J \right\} \quad (15)$$

that is, we are interested in solving the optimisation problem

$$\inf_{\lambda \in S_J, \Theta \in \mathbb{T}^J} \Psi_{\alpha}(\mu_{\lambda, \Theta} k),$$

with $J > 1$. Let us next denote $\lambda_n = (\lambda_{j,n})_{1 \le j \le J}$ and $\Theta_n = (\theta_{j,n})_{1 \le j \le J}$ for all $n \ge 1$. For convenience, we also introduce the shorthand notation $\mu_n k = \mu_{\lambda_n, \Theta_n} k$ and

$$\varphi_{j,n}^{(\alpha)}(y) = k(\theta_{j,n}, y) \left( \frac{\mu_n k(y)}{p(y)} \right)^{\alpha-1} \quad (16)$$

$$\dot{\varphi}_{j,n}^{(\alpha)} = \frac{\varphi_{j,n}^{(\alpha)}}{\int \varphi_{j,n}^{(\alpha)} d\nu}$$

for all $\alpha \in [0, 1)$, all $j = 1 \dots J$, all $n \ge 1$ and all $y \in \mathbb{Y}$.

Our goal in this section is to derive iterative schemes for the mixture weights and the mixture components parameters $(\lambda_n, \Theta_n)_{n \ge 1}$ ensuring that the sequence $(\Psi_{\alpha}(\mu_n k))_{n \ge 1}$ is decreasing. As Theorem 1 holds for any choice of parametric family, a first idea is to apply Theorem 1 to the variational family (15), which gives Corollary 3 below.

**Corollary 3** Assume (A1). Let $J \in \mathbb{N}^*$, let $\alpha \in [0, 1)$ and starting from an initial parameter set $(\lambda_1, \Theta_1) \in S_J^+ \times \mathbb{T}^J$, let $(\lambda_n, \Theta_n)_{n \ge 1}$ be defined iteratively such that for all $n \ge 1$,

$$\int_{\mathbb{Y}} \frac{(\mu_n k(y))^{\alpha} p(y)^{1-\alpha}}{\alpha - 1} \log \left( \frac{\mu_{n+1} k(y)}{\mu_n k(y)} \right) \nu(dy) \le 0. \quad (17)$$

Further assume that $\Psi_{\alpha}(\mu_1 k) < \infty$. Then, $\Psi_{\alpha}(\mu_{n+1} k) \le \Psi_{\alpha}(\mu_n k)$ for all $n \ge 1$.

Observe that we are in a less favourable situation in Corollary 3 with $J > 1$ compared to the cases we previously studied in Section 3. Indeed, we now have a ratio of sums inside the log function in (17), meaning that the approach from Theorem 1 to derive simple iterative schemes does not directly transfer to the variational family (15) for choices of kernel densities $k$ identified in Examples 1-2 and Lemma 1. However, by carefully exploiting the condition (17), we are able to overcome this difficulty in our second main theorem below.

**Theorem 2** Assume (A1). Let $J \in \mathbb{N}^*$, $\alpha \in [0, 1)$ and starting from an initial parameter set $(\lambda_1, \Theta_1) \in S_J^+ \times \mathbb{T}^J$, let $(\lambda_n, \Theta_n)_{n \ge 1}$ be defined iteratively such that for all $n \ge 1$,

$$\int_{\mathbb{Y}} \sum_{j=1}^{J} \lambda_{j,n} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \log \left( \frac{\lambda_{j,n+1}}{\lambda_{j,n}} \right) \nu(dy) \le 0 \quad (18)$$

$$\int_{\mathbb{Y}} \sum_{j=1}^{J} \lambda_{j,n} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \log \left( \frac{k(\theta_{j,n+1}, y)}{k(\theta_{j,n}, y)} \right) \nu(dy) \le 0. \quad (19)$$

Further assume that $\Psi_{\alpha}(\mu_1 k) < \infty$. Then, $\Psi_{\alpha}(\mu_{n+1} k) \le \Psi_{\alpha}(\mu_n k)$ for all $n \ge 1$.
**Proof** By Corollary 3, we can conclude if we show that (18) and (19) together imply (17). To show this, first observe that since the function $u \mapsto \frac{1}{\alpha-1} \log(u)$ is convex and $\alpha \in [0, 1)$, Jensen's inequality implies that: for all $y \in Y$ and all $n \ge 1$,

$$ \frac{1}{\alpha-1} \log \left( \frac{\mu_{n+1}k(y)}{\mu_n k(y)} \right) = \frac{1}{\alpha-1} \log \left( \sum_{j=1}^{J} \frac{\lambda_{j,n}k(\theta_{j,n}, y)}{\sum_{\ell=1}^{J} \lambda_{\ell,n}k(\theta_{\ell,n}, y)} \frac{\lambda_{j,n+1}k(\theta_{j,n+1}, y)}{\lambda_{j,n}k(\theta_{j,n}, y)} \right) \\ \le \sum_{j=1}^{J} \frac{\lambda_{j,n}k(\theta_{j,n}, y)}{\sum_{\ell=1}^{J} \lambda_{\ell,n}k(\theta_{\ell,n}, y)} \frac{1}{\alpha-1} \log \left( \frac{\lambda_{j,n+1}k(\theta_{j,n+1}, y)}{\lambda_{j,n}k(\theta_{j,n}, y)} \right) $$

that is:

$$ \frac{1}{\alpha-1} \log \left( \frac{\mu_{n+1}k(y)}{\mu_n k(y)} \right) \le \sum_{j=1}^{J} \frac{\lambda_{j,n}}{\alpha-1} \frac{k(\theta_{j,n}, y)}{\mu_n k(y)} \log \left( \frac{\lambda_{j,n+1}k(\theta_{j,n+1}, y)}{\lambda_{j,n}k(\theta_{j,n}, y)} \right). $$

Multiplying by $(\mu_n k(y))^\alpha p(y)^{1-\alpha}$ on both sides, integrating with respect to $\nu(dy)$ and using the definition of $\varphi_{j,n}^{(\alpha)}(y)$ in (16), this in turn implies that: for all $n \ge 1$,

$$ \int_Y \frac{(\mu_n k(y))^\alpha p(y)^{1-\alpha}}{\alpha-1} \log \left( \frac{\mu_{n+1}k(y)}{\mu_n k(y)} \right) \nu(dy) \\ \le \int_Y \sum_{j=1}^J \lambda_{j,n} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \log \left( \frac{\lambda_{j,n+1}k(\theta_{j,n+1}, y)}{\lambda_{j,n}k(\theta_{j,n}, y)} \right) \nu(dy). $$

As a consequence, the condition

$$ \int_Y \sum_{j=1}^{J} \lambda_{j,n} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \log \left( \frac{\lambda_{j,n+1}k(\theta_{j,n+1}, y)}{\lambda_{j,n}k(\theta_{j,n}, y)} \right) \nu(dy) \le 0 \quad (20) $$

implies (17). The condition (20) in itself is then straightforwardly implied by (18) and (19) and the proof is concluded. ■

Strikingly, (18) does not depend on $\Theta_{n+1}$ nor does (19) depend on $\lambda_{n+1}$ in Theorem 2. This means that we can treat these two conditions separately and thus that the weights and components parameters of the mixture can be optimised simultaneously. This result was far from immediate by looking at the initial condition (17) and, as we shall see later in Section 8, will lead to reduced computational power in practice.

In addition, we have recovered in (19) the key property used in Section 3: compared to the condition (17) which involved the log of weighted sums of kernel densities, (19) considers a sum of logs of kernel densities. This suggests that we can extend the updates derived in Section 3 for well-chosen kernel densities $k$ to the more general case of mixture models.

Observe finally that the dependency in $\lambda_{j,n+1}$ appearing in (18) is simpler than the dependency in $\theta_{j,n+1}$ appearing in (19) and that is expressed through the kernel density $k$. As a result, we will first study (18), in the hope of deriving iterative update formulas for the mixture weights that do not require a specific choice of kernel density $k$. Interestingly, while the natural idea is to perform direct optimisation of the left-hand side of (18), we will derive a more general expression for the mixture weights update, which shall induce numerical advantages later illustrated in Section 8.
## 4.1 Choice of $(\lambda_n)_{n \ge 1}$

In the following theorem, we identify an update formula which satisfies (18), regardless of the choice of the kernel density $k$.

**Theorem 3** Assume (A1). Let $\alpha \in [0, 1)$, let $(\eta_n)_{n \ge 1}$ be valued in $(0, 1]$, let $(\kappa_n)_{n \ge 1}$ be such that $(\alpha - 1)\kappa_n \ge 0$ at all times $n \ge 1$ and let $(\Theta_n)_{n \ge 1}$ be any sequence valued in $\mathbb{T}^J$. Starting from an initial $\lambda_1 \in S_J^+$, let $(\lambda_n)_{n \ge 1}$ be defined iteratively such that for all $n \ge 1$

$$ \lambda_{j,n+1} = \frac{\lambda_{j,n} \left[ \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1)\kappa_n \right]^{\eta_n}}{\sum_{\ell=1}^{J} \lambda_{\ell,n} \left[ \int_Y \varphi_{\ell,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1)\kappa_n \right]^{\eta_n}}, \quad j = 1 \dots J. \quad (21) $$

Then (18) holds.

**Proof** We first check that the integrals appearing in (21) are finite:

$$ \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy) < \infty, \quad j = 1, \dots, J, \quad n \ge 1. \quad (22) $$

Using (16), that $\lambda_{j,n} > 0$, that $\lambda_{j,n}k(\theta_{j,n}, y) \le \mu_n k(y)$ and Jensen's inequality applied to the concave function $u \mapsto u^{1-\alpha}$, we have, for any $j = 1, \dots, J$ and $n \ge 1$,

$$ \begin{aligned} \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy) &= \int_Y k(\theta_{j,n}, y) \left( \frac{\mu_n k(y)}{p(y)} \right)^{\alpha-1} \nu(dy) \\ &\le \frac{1}{\lambda_{j,n}} \int_Y \mu_n k(y) \left( \frac{p(y)}{\mu_n k(y)} \right)^{1-\alpha} \nu(dy) \\ &\le \frac{1}{\lambda_{j,n}} \left( \int_Y p(y) \nu(dy) \right)^{1-\alpha}. \end{aligned} $$

The bound (22) follows by (A1). Now, to prove (18), we treat the cases $\eta_n = 1$ and $\eta_n \in (0, 1)$ separately.

(a) Case $\eta_n = 1$. Since $(\alpha - 1)\kappa_n \ge 0$ with $\alpha \in (0, 1)$, we have that

$$ \kappa_n \sum_{j=1}^{J} \lambda_{j,n} \log(\lambda_j / \lambda_{j,n}) \ge 0 $$

where we have used that $\sum_{j=1}^{J} \lambda_{j,n} \log(\lambda_j / \lambda_{j,n}) \le \sum_{j=1}^{J} \lambda_{j,n}(\lambda_j / \lambda_{j,n} - 1) = 0$. In other words, to obtain (18) in the particular case $\eta_n = 1$, it is enough to show

$$ \int_Y \sum_{j=1}^{J} \lambda_{j,n} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \log \left( \frac{\lambda_{j,n+1}}{\lambda_{j,n}} \right) \nu(dy) + \kappa_n \sum_{j=1}^{J} \lambda_{j,n} \log \left( \frac{\lambda_{j,n+1}}{\lambda_{j,n}} \right) \le 0 $$

that is, since (22) holds,

$$ \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \nu(dy) + \kappa_n \right] \log \left( \frac{\lambda_{j,n+1}}{\lambda_{j,n}} \right) \le 0. \quad (23) $$
Notice then that by definition of $(\lambda_{j,n+1})_{1\le j\le J}$ when $\eta_n = 1$, we can write

$$ \lambda_{n+1} = \underset{\lambda \in S_j^+}{\operatorname{argmin}} \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \log \left( \frac{\lambda_j}{\lambda_{j,n}} \right). $$

[Indeed, setting $\beta_j = \lambda_{j,n} \left[ \int_Y \varphi_{j,n}^{(\alpha)}(y)\nu(dy) + (\alpha-1)\kappa_n \right]$ and $\bar{\beta}_j = \beta_j / \sum_{\ell=1}^J \beta_\ell$ for all $j = 1...J$, we have that $\sum_{j=1}^J \bar{\beta}_j \log(\bar{\beta}_j/\lambda_j) \ge 0$ and that this quantity is minimal when $\lambda_j = \bar{\beta}_j$ for $j = 1...J$.] This implies (23) and settles the case $\eta_n = 1$.

(b) For the particular case $\eta_n \in (0,1)$, we will use that for all $\epsilon > 0$ and all $u > 0$,

$$ \log(u) = \frac{1}{\epsilon} \log(u^{\epsilon}) \ge \frac{1}{\epsilon} \left(1 - \frac{1}{u^{\epsilon}}\right). $$

Indeed, since $\int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \le 0$ for all $j = 1...J$, we can then write that for all $\epsilon > 0$,

$$ \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_{Y} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \log \left( \frac{\lambda_{j,n+1}}{\lambda_{j,n}} \right) \\ \le \frac{1}{\epsilon} \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_{Y} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \left[ 1 - \left( \frac{\lambda_{j,n}}{\lambda_{j,n+1}} \right)^{\epsilon} \right]. \quad (24) $$

Now notice that by definition of $(\lambda_{j,n+1})_{1\le j\le J}$ we can write

$$ \lambda_{n+1} = \underset{\lambda \in S_j^+}{\operatorname{argmin}} \frac{1}{\epsilon} \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \left[ 1 - \left( \frac{\lambda_{j,n}}{\lambda_j} \right)^{\epsilon} \right] $$

when $\epsilon$ satisfies $\eta_n = \frac{1}{1+\epsilon}$. [Indeed setting $\beta_j = \lambda_{j,n} \left[ \int_Y \varphi_{j,n}^{(\alpha)}(y)\nu(dy) + (\alpha-1)\kappa_n \right]^{\frac{1}{1+\epsilon}}$ and $\bar{\beta}_j = \beta_j / \sum_{\ell=1}^J \beta_\ell$ for all $j = 1...J$, we have by convexity of the function $u \mapsto u^{1+\epsilon}$ that $\sum_{j=1}^J (\bar{\beta}_j/\lambda_j)^{1+\epsilon} \lambda_j \ge (\sum_{j=1}^J \bar{\beta}_j)^{1+\epsilon}$ and that this quantity is minimal when $\lambda_j = \bar{\beta}_j$ for $j = 1...J$.] We then deduce that taking $\epsilon = \eta_n^{-1} - 1$ (it is always possible since $\eta_n \in (0,1)$ by assumption) yields

$$ \frac{1}{\epsilon} \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_{Y} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \left[ 1 - \left( \frac{\lambda_{j,n}}{\lambda_{j,n+1}} \right)^{\epsilon} \right] \le 0 $$

which in turn yields (18) [since combined with (24) it implies (23) which itself implies (18) as seen in the case $\eta_n = 1$. This settles the case $\eta_n \in (0,1)$.

<figure><img src="image_4.png" /></figure>

Notice that as a byproduct of the proof of Theorem 3, the mixture weights update given by (21) can be rewritten under the form: for all $n \ge 1$,

$$ \lambda_{n+1} = \underset{\lambda \in S_j^+}{\operatorname{argmin}} h_n(\lambda) $$
where, setting $\epsilon = \eta_n^{-1} - 1$, we have defined for all $\lambda \in S_J^+$,

$$h_n(\lambda) = \begin{cases} \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \log \left( \frac{\lambda_j}{\lambda_{j,n}} \right), & \text{if } \eta_n = 1, \\ \frac{1}{\epsilon} \sum_{j=1}^{J} \lambda_{j,n} \left[ \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \nu(dy) + \kappa_n \right] \left[ 1 - \left( \frac{\lambda_{j,n}}{\lambda_j} \right)^{\epsilon} \right], & \text{if } \eta_n \in (0, 1). \end{cases}$$

More specifically, $h_n(\lambda)$ acts as an upper bound of the left-hand side of (18) and we recover exactly the left-hand side of (18) in the particular case $\eta_n = 1$ and $\kappa_n = 0$.

Now that we have established the mixture weights updates (21) in Theorem 3, we are interested in deriving update formulas for the sequence $(\Theta_n)_{n \ge 1}$ satisfying (19), which we will then pair up with (21) in order to apply Theorem 2. From this point onwards, all the proofs of the coming results will be deferred to the appendices to ease the reading.

## 4.2 Choice of $(\Theta_n)_{n \ge 1}$

We investigate two different approaches for choosing $(\Theta_n)_{n \ge 1}$.

### 4.2.1 A MAXIMISATION APPROACH

As done in Corollary 1, an idea is to consider an update for $(\Theta_n)_{n \ge 1}$ of the form

$$\Theta_{n+1} = \underset{\Theta \in \mathcal{T}^J}{\operatorname{argmax}} g_n(\Theta), \quad n \ge 1,$$

where the function $g_n$ is constructed as a lower bound on $\Theta \in \mathcal{T}^J$ of the function $\Theta \mapsto \int_Y \sum_{j=1}^J \lambda_{j,n} \varphi_{j,n}^{(\alpha)}(y) \log(k(\theta_j, y)/k(\theta_{j,n}, y)) \nu(dy)$ that satisfies $g_n(\Theta_n) = 0$. In doing so, the function $\tilde{g}_n : \Theta \mapsto g_n(\Theta)/(\alpha - 1)$ evaluated at $\Theta_{n+1}$ is an upper bound of the left-hand side of (19) and $\tilde{g}_n(\Theta_{n+1}) \le 0$ implies (19). This leads us to Corollary 4 below.

**Corollary 4 (Generalised maximisation approach)** Assume (A1). Let $\alpha \in [0, 1)$, let $(\eta_n)_{n \ge 1}$ be valued in $(0, 1]$ and let $(\kappa_n)_{n \ge 1}$ be such that $(\alpha - 1)\kappa_n \ge 0$ at all times $n \ge 1$. Furthermore, let $(b_{j,n})_{n \ge 1}$ be a non-negative sequence for all $j = 1...J$. Starting from an initial parameter set $(\lambda_1, \Theta_1) \in S_J^+ \times \mathcal{T}^J$, let $(\lambda_n, \Theta_n)_{n \ge 1}$ be defined iteratively for all $n \ge 1$ in such a way that (21) holds and

$$\theta_{j,n+1} = \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_Y \left[ \varphi_{j,n}^{(\alpha)}(y) + b_{j,n}k(\theta_{j,n}, y) \right] \log \left( \frac{k(\theta, y)}{k(\theta_{j,n}, y)} \right) \nu(dy), \quad j = 1...J, \quad (25)$$

where we assume that this argmax is uniquely defined at each step. Then, we can apply Theorem 2.

The proof of this result is deferred to Appendix B.1. Under the assumptions of Corollary 4, Algorithm 1 leads to a systematic decrease in $\Psi_\alpha$ at each step. This result effectively generalises the monotonicity property from Corollary 1 to the case of mixture models and we can deduce simple iterative schemes satisfying (25) for well-chosen kernel densities $k$. To illustrate this, we provide below the extension of Example 1 to Gaussian Mixture Models (and Lemma 1 can be extended in the same way, see Appendix B.2).
**Algorithm 1:** Maximisation approach algorithm for mixture models

At iteration n,

For all j = 1...J, set

$$
\lambda_{j,n+1} = \frac{\lambda_{j,n} \left[ \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1) \kappa_n \right]^{\eta_n}}{\sum_{\ell=1}^{J} \lambda_{\ell,n} \left[ \int_Y \varphi_{\ell,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1) \kappa_n \right]^{\eta_n}}
$$
$$
\theta_{j,n+1} = \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_Y \left[ \varphi_{j,n}^{(\alpha)}(y) + b_{j,n} k(\theta_{j,n}, y) \right] \log \left( \frac{k(\theta, y)}{k(\theta_{j,n}, y)} \right) \nu(dy).
$$

**Example 3 (Maximisation approach for Gaussian Mixture Models)** We consider the case of d-dimensional Gaussian mixture densities, in which case $Y = \mathbb{R}^d$, $\nu$ is the d-dimensional Lebesgue measure and $k(\theta_j, y) = \mathcal{N}(y; m_j, \Sigma_j)$, where $\theta_j = (m_j, \Sigma_j) \in \mathcal{T}$ denotes the mean and covariance matrix of the j-th Gaussian component density. We let the parameter set $\mathcal{T}$ include all possible means in $\mathbb{R}^d$ and all possible positive-definite covariance matrices. Finally, we assume that the non-negative function $p$ defined on $Y$ satisfies (8).

Then (A1) holds. Moreover, starting from any $(\lambda_1, \Theta_1) \in S_J^+ \times \mathcal{T}^J$ so that $\Psi_\alpha(\mu_1 k) < \infty$
and denoting $\theta_{j,n} = (m_{j,n}, \Sigma_{j,n})$ for all $j = 1 \dots J$ and all $n \ge 1$, it holds that: for all
$j = 1 \dots J$, all $n \ge 1$ and all $\gamma_{j,n} \in (0,1]$, there exists $b_{j,n} \ge 0$ such that the argmax problem
(25) has a unique solution defined by:

$$
m_{j,n+1} = \gamma_{j,n} \hat{m}_{j,n} + (1 - \gamma_{j,n}) m_{j,n} \\
\Sigma_{j,n+1} = \gamma_{j,n} \hat{\Sigma}_{j,n} + (1 - \gamma_{j,n}) \Sigma_{j,n} + \gamma_{j,n} (1 - \gamma_{j,n}) (\hat{m}_{j,n} - m_{j,n}) (\hat{m}_{j,n} - m_{j,n})^T,
$$

where, using the definition of $\check{\varphi}_{j,n}^{(\alpha)}$ in (16), we set

$$
\begin{align*}
\hat{m}_{j,n} &= \int_Y y \check{\varphi}_{j,n}^{(\alpha)}(y) \nu(dy), \\
\hat{\Sigma}_{j,n} &= \int_Y y y^T \check{\varphi}_{j,n}^{(\alpha)}(y) \nu(dy) - \hat{m}_{j,n} \hat{m}_{j,n}^T.
\end{align*}
$$

The detailed derivation of the mean and covariance updates above is deferred to Section 6.2.
The interpretation of these updates is akin to the one already made in Example 1: $b_{j,n}$ acts
as a regularisation parameter which, through $\gamma_{j,n}$, permits a tradeoff between an update close
to the current parameter $\theta_{j,n}$ and choosing the Gaussian density with exactly the same mean
and covariance matrix as $\check{\varphi}_{j,n}^{(\alpha)}$. As per written in Corollary 4, these updates are compatible
with the mixture weights updates (21), resulting in a systematic decrease of $(\Psi_{\alpha}(\mu_n k))_{n \ge 1}$.

We next present another possible update formula for $(\lambda_n, \Theta_n)_{n \ge 1}$.

4.2.2 A GRADIENT-BASED APPROACH

In the spirit of Corollary 2, we now resort to Gradient Descent steps to satisfy (18).
**Algorithm 2:** Gradient-based approach algorithm for mixture models

At iteration $n$,

For all $j = 1...J$, set

$$ 
\begin{aligned}
\lambda_{j,n+1} &= \frac{\lambda_{j,n} \left[ \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1)\kappa_n \right]^{\eta_n}}{\sum_{\ell=1}^{J} \lambda_{\ell,n} \left[ \int_Y \varphi_{\ell,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1)\kappa_n \right]^{\eta_n}} \\
\theta_{j,n+1} &= \theta_{j,n} - \frac{\gamma_{j,n}}{\beta_{j,n}} \nabla g_{j,n}(\theta_{j,n}).
\end{aligned}
$$

**Corollary 5 (Generalised gradient-based approach)** Assume (A1). Let $\mathcal{T} \subseteq \mathbb{R}^d$ be a convex set, let $\alpha \in [0, 1]$, let $(\eta_n)_{n \ge 1}$ be valued in $(0, 1]$ and let $(\kappa_n)_{n \ge 1}$ be such that $(\alpha - 1)\kappa_n \ge 0$ at all times $n$. Furthermore, for all $j = 1...J$, let $(\gamma_{j,n})_{n \ge 1}$ be valued in $(0, 1]$. Starting from an initial parameter set $(\lambda_1, \Theta_1) \in S_J^+ \times \mathcal{T}^J$, let $(\lambda_n, \Theta_n)_{n \ge 1}$ be defined iteratively for all $n \ge 1$ in such a way that (21) holds and

$$ \theta_{j,n+1} = \theta_{j,n} - \frac{\gamma_{j,n}}{\beta_{j,n}} \nabla g_{j,n}(\theta_{j,n}), \quad j = 1...J, \qquad (26) $$

where for all $j = 1...J$, $(g_{j,n})_{n \ge 1}$ is defined by: for all $n \ge 1$ and all $\theta \in \mathcal{T}$,

$$ g_{j,n}(\theta) = \int_{\mathcal{Y}} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \log \left( \frac{k(\theta, y)}{k(\theta_{j,n}, y)} \right) \nu(dy). $$

and $g_{j,n}$ is assumed to be $\beta_{j,n}$-smooth. Then, we can apply Theorem 2.

The proof of this result is deferred to Appendix B.3. Under the assumptions of Corollary 5, Algorithm 2 ensures a systematic decrease in $\Psi_\alpha$ at each step and Corollary 5 thus extends Corollary 2 to mixture models. Much like what we did for Corollary 2, we want to identify how our updates relate to Gradient Descent-based techniques for optimising $\Theta$. Under common differentiability assumptions, we have: for all $n \ge 1$, all $j = 1...J$ and all $\theta' \in \mathcal{T}$,

$$ \nabla g_{j,n}(\theta') = \int_{\mathcal{Y}} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta,y)=(\theta',y)} \nu(dy), $$

so that (26) becomes

$$ \theta_{j,n+1} = \theta_{j,n} - \frac{\gamma_{j,n}}{\beta_{j,n}} \int_{\mathcal{Y}} \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta,y)=(\theta_{j,n},y)} \nu(dy), \quad j = 1 \dots J. \qquad (27) $$

The link with Gradient Descent-based Variational Inference shall become apparent by (i) writing the update formulas that ensue from Gradient Descent iterations for $\alpha$-divergence and Rényi's $\alpha$-divergence minimisation and (ii) understanding how Example 2 generalises to Gaussian Mixture Models. Given $(\lambda_n, \Theta_n) \in S_J^+ \times \mathcal{T}^J$, an index $j$ in $1...J$ and letting
p = p(·, D), performing one Gradient Descent iteration w.r.t. θ<sub>j,n</sub> for α-divergence and Rényi's α-divergence minimisation indeed respectively amounts to updating θ<sub>j,n+1</sub> as follows

$$
\begin{align*}
\theta_{j,n+1} &= \theta_{j,n} - r_{j,n} \lambda_{j,n} \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta,y)=(\theta_{j,n},y)} \nu(dy), \\
\theta_{j,n+1} &= \theta_{j,n} - r_{j,n} \frac{\lambda_{j,n}}{\sum_{j=1}^J \lambda_{j,n} \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy)} \int_Y \frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta,y)=(\theta_{j,n},y)} \nu(dy)
\end{align*}
$$

where $r_{j,n} > 0$ is the learning rate (we refer to Appendix B.4 for details regarding these
updates). There is thus a similarity between (27) and the Gradient Descent updates above.
To fully comprehend the connection between our gradient-based approach and Gradient
Descent steps, we present below the generalisation of Example 2 to Gaussian Mixture Models.
The smoothness assumption on $g_{j,n}$ will be satisfied in that example and we will recover a
Gradient Descent scheme for Rényi's $\alpha$-divergence minimisation as a special case.

**Example 4 (Gradient-based approach for Gaussian Mixture Models)** We consider the case of d-dimensional Gaussian mixture densities with $k(\theta_j, y) = \mathcal{N}(y; m_j, \Sigma_j)$, where this time $\theta_j = m_j$. Further assume that $\mathcal{T} \subseteq \mathbb{R}^d$ is a convex subset and that $\Sigma_j = \sigma_j^2 \mathbf{I}_d$, where $\sigma_j^2 > 0$ and $\mathbf{I}_d$ denotes the d-dimensional identity matrix. Finally, we assume that the non-negative function $p$ defined on $\mathcal{Y}$ satisfies (13).

Then, (A1) holds. Setting $\beta_{j,n} = \sigma_j^{-2}(1-\alpha)^{-1} \int_Y \varphi_{j,n}^{(\alpha)}(y)\nu(dy)$ and denoting $\theta_{j,n} = m_{j,n}$, the function $g_{j,n}$ is $\beta_{j,n}$-smooth for all $j = 1 \dots J$ and all $n \ge 1$ so that the means of the Gaussian densities $(\mathcal{N}(y; m_j, \Sigma_j))_{1 \le j \le J}$ can be optimised as follows: for all $n \ge 1$,

$$
m_{j,n+1} = \gamma_{j,n} \hat{m}_{j,n} + (1 - \gamma_{j,n}) m_{j,n}, \quad j = 1 \dots J,
$$

where $\check{\varphi}_{j,n}^{(\alpha)}$ is defined in (16), $\gamma_{j,n} \in [0,1)$ and $\hat{m}_{j,n} = \int_Y y \check{\varphi}_{j,n}^{(\alpha)}(y) \nu(dy)$ (we refer to Corollary 8 for details). In particular, letting $(\gamma'_{j,n})_{n \ge 1}$ be valued in $(0,1]$ for all $j=1 \dots J$,

$$
\gamma_{j,n} = \gamma'_{j,n} \frac{\lambda_{j,n} \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy)}{\int_Y \mu_n k(y)^{\alpha} p(y)^{1-\alpha} \nu(dy)} \in (0,1], \quad j=1 \dots J, \quad n \ge 1,
$$

(since $\int_Y \mu_n k(y)^\alpha p(y)^{1-\alpha} \nu(dy) = \sum_{j=1}^J \lambda_{j,n} \int_Y \varphi_{j,n}^{(\alpha)}(y) \nu(dy)$). The resulting iterative algorithm is given by the following update at time $n \ge 1$:

$$
m_{j,n+1} = m_{j,n} + \gamma'_{j,n} \frac{\int_Y \lambda_{j,n} \varphi_{j,n}^{(\alpha)}(y) (y - m_{j,n}) \nu(dy)}{\int_Y \mu_n k(y)^\alpha p(y)^{1-\alpha} \nu(dy)}, \quad j = 1 \dots J. \quad (28)
$$

Interestingly, (28) can be seen as a Gradient Descent step w.r.t. $\theta_{j,n}$ for Rényi's $\alpha$-divergence minimisation with a learning rate $r_{j,n} = \sigma_j^2(1-\alpha)\gamma'_{j,n}$. Hence, if we were to solely rely on the Gradient Descent literature, the convergence of the iterative sequence $(\theta_{j,n})_{n\ge 1}$ defined by (28) would require the sequence $(\lambda_n)_{n\ge 1}$ to be constant. This is in contrast with Corollary 5, which allows for a simultaneous optimisation of $\lambda$ and $\Theta$ according to (21) and (28).

We now add on the comments made in Section 3 for the gradient-based methodology (which
we built in Corollary 2 and have since extended to mixture models in Corollary 5):
* A core insight from Corollary 5, which is exemplified in Example 4, is that under a smoothness assumption on $g_{j,n}$ our mixture weights iterative updates are compatible with gradient-based updates, themselves linked to the Gradient Descent literature. In other words, we have embedded Gradient Descent-based iterative updates, which only act on $\Theta$, within a larger framework where simultaneous updates for the mixture components parameters and the mixture weights are well-supported theoretically.

* Putting things into perspective with the maximisation approach once again, notice that our previous conclusions from Section 3 still hold. Namely: (i) the updates on the means from Examples 3 and 4 coincide, meaning that contrary to the gradient-based approach, the maximisation approach enables covariance matrices updates on top of means updates (ii) the maximisation approach permits us to bypass the smoothness assumption made in Corollary 5. As a whole, these properties make the maximisation approach a compelling alternative to the gradient-based approach.

We presented two approaches to construct iterative schemes $(\lambda_n, \Theta_n)_{n \ge 1}$ that ensure a monotonic decrease in $\Psi_\alpha$ at each step and lead to simple updates formulas for Gaussian Mixture Models. We now describe how our framework can be linked to the existing literature.

# 5. Related Work

In this section, we detail how our work relates to and improve on previous algorithms proposed for $\alpha$-/Rényi's $\alpha$-divergence minimisation.

## 5.1 Rényi Divergence Variational Inference (Li and Turner, 2016)

In Li and Turner (2016), they seek to maximise the Variational Rényi (VR) Bound via (Stochastic) Gradient Ascent. This objective function is derived from Rényi's $\alpha$-divergence and is given by: for all variational density $q \in \mathcal{Q}$ and all $\alpha \in \mathbb{R} \setminus \{1\}$,

$$ \mathcal{L}_{\alpha}(q, \mathcal{D}) := \frac{1}{1-\alpha} \log \left( \int_{\mathbb{Y}} q(y)^{\alpha} p(y, \mathcal{D})^{1-\alpha} \nu(dy) \right). \quad (29) $$

Contrary to us, the work from Li and Turner (2016) does not consider the case where $q$ belongs to the mixture models variational family (15), that is $q = \mu_{\lambda, \Theta} k$.

Yet, a parallel can be drawn in the GMM case between Li and Turner (2016) and our approach by observing that the gradient-based updates on the means in Example 4 each coincide with a Gradient Ascent step on the objective function (29) for a well-chosen learning rate (since the gradient of $\mathcal{L}_\alpha$ is proportional to the gradient of Rényi's $\alpha$-divergence with a factor $-\alpha^{-1}$, this follows from the remarks made in Section 4.2.2 regarding Gradient Descent steps for Rényi's $\alpha$-divergence minimisation in the GMM case).

Our work hence provides a theoretical framework which enables simultaneous optimisation of the mixture weights $\lambda$ and of the mixture components parameters $\Theta$. In addition, beyond the gradient-based updates, we propose the novel maximisation updates and we allow for covariance matrices optimisation. We also emphasise that our maximisation approach will apply to well-chosen kernels $k$ beyond the Gaussian case, as we will detail in Section 6.
## 5.2 The Power Descent Algorithm (Daudel et al., 2021)

In order to identify the connection between our work and the Power Descent algorithm introduced in Daudel et al. (2021), let us briefly present the latter.

The Power Descent is a gradient-based algorithm that operates on probability measures and performs α-divergence minimisation for all $\alpha \in \mathbb{R} \setminus \{1\}$. More precisely, equipping $\mathcal{T}$ with a $\sigma$-field $\mathcal{T}$ and denoting by $M_1(\mathcal{T})$ the space of probability measures on $(\mathcal{T}, \mathcal{T})$, the Power Descent optimises $\Psi_\alpha(\mu k)$ with respect to $\mu \in M_1(\mathcal{T})$, where $\mu k(y) = \int_{\mathcal{T}} \mu(d\theta)k(\theta, y)$ for all $\mu \in M_1(\mathcal{T})$ and all $y \in \mathcal{Y}$. Given an initial probability measure $\mu_1 \in M_1(\mathcal{T})$, it does so by performing several one-step transitions of the Power Descent algorithm:

$$ \mu_{n+1} = \mathcal{I}_{\alpha}(\mu_n), \quad n \ge 1, \qquad (30) $$

where, for all $\mu \in M_1(\mathcal{T})$ and all $\theta \in \mathcal{T}$,

$$ b_{\mu,\alpha}(\theta) = \int_{\mathcal{Y}} k(\theta, y) \frac{1}{\alpha - 1} \left[ \left( \frac{\mu k(y)}{p(y)} \right)^{\alpha - 1} - 1 \right] \nu(dy) $$

$$ \mathcal{I}_{\alpha}(\mu)(d\theta) = \frac{\mu(d\theta) [(\alpha - 1)(b_{\mu,\alpha}(\theta) + \kappa) + 1]^{\frac{\eta}{1-\alpha}}}{\mu([\alpha - 1](b_{\mu,\alpha}(\theta) + \kappa) + 1]^{\frac{\eta}{1-\alpha}})} $$

Daudel et al. (2021) motivated the Power Descent algorithm by establishing a monotonicity result for this algorithm obtained as a particular case of (Daudel et al., 2021, Theorem 1). In their result, the monotonic decrease in $\Psi_\alpha$ of the scheme (30) holds for all $\mu \in M_1(\mathcal{T})$, all $\alpha \in \mathbb{R} \setminus \{1\}$, all $\eta \in (0, 1]$ and all $\kappa$ such that $(\alpha - 1)\kappa \ge 0$. We provide below a more general version of their result, where the monotonic decrease in $\Psi_\alpha$ holds for well-chosen values of $\eta$ that are larger than 1 when $\alpha < 0$ (we refer to Appendix B.5 for the proof of this result).

**Proposition 2** Assume that $p$ and $k$ are as in (A1). Let $(\alpha, \eta)$ belong to any of the following cases.

(i) $\alpha \le -1$ and $\eta \in (0, (\alpha - 1)/\alpha]$;

(ii) $\alpha \in (-1, 0)$ and $\eta \in (0, 1 - \alpha]$;

(iii) $\alpha \in [0, 1)$ or $\alpha > 1$ and $\eta \in (0, 1]$.

Moreover, let $\mu \in M_1(\mathcal{T})$ be such that $\Psi_\alpha(\mu k) < \infty$ and let $\kappa$ be such that $(\alpha - 1)\kappa \ge 0$. Then, the two following assertions hold.

(i) We have $\Psi_\alpha(\mathcal{I}_\alpha(\mu)k) \le \Psi_\alpha(\mu k)$.

(ii) We have $\Psi_\alpha(\mathcal{I}_\alpha(\mu)k) = \Psi_\alpha(\mu k)$ if and only if $\mu = \mathcal{I}_\alpha(\mu)$.

Building on the monotonicity result provided by (Daudel et al., 2021, Theorem 1) for the Power Descent algorithm - that we generalised in Proposition 2 - Daudel et al. (2021) then applied this algorithm to mixture weights optimisation by letting the initial probability measure $\mu_1 \in M_1(\mathcal{T})$ be a weighted sum of Dirac measures of the form $\mu_1 = \sum_{j=1}^J \lambda_{j,n} \delta_{\theta_j}$, with
$\Theta \in \mathcal{T}$ and $\lambda_1 \in S_J$. For that choice of $\mu_1, \mu_n$ in (30) can be written as $\mu_n = \sum_{j=1}^J \lambda_{j,n} \delta_{\theta_j}$ at time $n$ and the Power Descent amounts to performing the update

$$ \lambda_{j,n+1} = \frac{\lambda_{j,n} \left[ (\alpha - 1)(b_{\mu_n,\alpha}(\theta_j) + \kappa) + 1 \right]^{\frac{\eta}{1-\alpha}}}{\sum_{\ell=1}^{J} \lambda_{\ell,n} \left[ (\alpha - 1)(b_{\mu_n,\alpha}(\theta_\ell) + \kappa) + 1 \right]^{\frac{\eta}{1-\alpha}}}, \quad j = 1 \dots J. \quad (31) $$

Interestingly, the update (31) corresponds to the update on the mixture weights (21) we have identified in Theorem 3 for $\Theta_n = \Theta$, $\eta_n = \eta/(1-\alpha)$ and $\kappa_n = \kappa$. Steaming from this link between our approach and Daudel et al. (2021), we now make two important comments.

* *Benefits of our approach compared to Daudel et al. (2021).* The monotonicity result from Daudel et al. (2021) generalised in Proposition 2 requires the sequence $(\Theta_n)_{n \ge 1}$ to be constant when applied to mixture weights optimisation in (31). This restricts the variational family to mixture models with fixed mixture components parameters. To remedy this problem, Daudel et al. (2021) proposed a fully-adaptive algorithm that alternates between an Exploitation step optimising the mixture weights according to (31) and an Exploration step acting on the mixture components parameters. However, they established no theoretical guarantees for their complete Exploitation-Exploration algorithm, as the choice of the Exploration step remained mostly unexplored.

A strong improvement of our approach is then that we provide theoretically-sound updates for $(\lambda_n, \Theta_n)_{n \ge 1}$, with the particularity that our mixture weights updates relate to the framework of Daudel et al. (2021). In that sense, our approach supplements the work done in Daudel et al. (2021) (albeit by an entirely different proof technique). Furthermore, we do not need to alternate between mixture weights and mixture components parameters updates in our algorithms, as both can be carried out simultaneously (as done in Algorithms 1 and 2). In practice, this will permit us to reduce the computational cost (as the samples will be shared throughout the mixture weights and the mixture components parameters approximated updates, see Section 7).

* *A gradient-based mixture weights update.* Daudel et al. (2021) establishes that the Power Descent belongs to a family of gradient-based algorithms which includes the Entropic Mirror Descent algorithm, a typical optimisation algorithm for optimisation under simplex constraints, as a special case. Viewed from this angle, the parameter $\eta$ in (31) can be understood as a learning rate with $b_{\mu_n,\alpha}$ playing the role of the gradient of $\Psi_\alpha$. Connecting the Power Descent to our framework thus sheds light on the gradient-based nature of the mixture weights update (21) from Theorem 3 and provides a better understanding of the role of the parameter $\eta_n$ appearing in this update. This aspect will be helpful to interpret our numerical experiments in Section 8.

## 5.3 The M-PMC Algorithm (Cappé et al., 2008)

For any measurable positive function $p$ on $(\mathcal{Y}, \mathcal{Y})$, the M-PMC algorithm (Cappé et al., 2008) aims at solving the optimisation problem

$$ \sup_{(\lambda \in S_J, \Theta \in \mathcal{T}^J)} \int_{\mathcal{Y}} \log \left( \sum_{j=1}^{J} \lambda_j k(\theta_j, y) \right) p(y) \nu(dy), $$
or equivalently, at minimising the inclusive Kullback-Leibler divergence $D_0(\mu_{\lambda,\Theta}K||\mathbb{P})$ w.r.t. to $(\lambda, \Theta)$, where $\mathbb{P}(A) = \int_A p \, d\nu / \int_Y p \, d\nu$ for all $A \in \mathcal{Y}$. This is done in (Cappé et al., 2008, Section 2) by introducing the following iterative updates: for all $j = 1...J$ and all $n \ge 1$,

$$ 
\begin{aligned}
\lambda_{j,n+1} &= \int_Y \frac{\lambda_{j,n} k(\theta_{j,n}, y)}{\sum_{\ell=1}^{J} \lambda_{\ell,n} k(\theta_{\ell,n}, y)} \frac{p(y)}{\int_Y p(y) \nu(dy)} \nu(dy) \\
\theta_{j,n+1} &= \underset{\theta_j \in \mathcal{T}}{\operatorname{argmax}} \int_Y \frac{\lambda_{j,n} k(\theta_{j,n}, y)}{\sum_{\ell=1}^{J} \lambda_{\ell,n} k(\theta_{\ell,n}, y)} \log(k(\theta_j, y)) p(y) \nu(dy).
\end{aligned}
$$

Cappé et al. (2008) motivated the two updates above by noticing that they can be seen as integrated versions under the target distribution of the update formulas for the Expectation-Maximisation (EM) algorithm applied to the mixture density parameter estimation problem

$$ \sup_{(\lambda \in \mathcal{S}_J, \Theta \in \mathcal{T}^J)} \sum_{m=1}^{M} \log \left( \sum_{j=1}^{J} \lambda_j k(\theta_j, Y_m) \right), $$

meaning that these updates ensure a systematic increase in the integrated likelihood at each step. As these updates also correspond to the case $\alpha = 0$, $\eta_n = 1$, $\kappa_n = 0$, $a_{j,n} = 1$ and $b_{j,n} = 0$ in Algorithm 1, the M-PMC algorithm is in fact included in our framework and we can interpret our theoretical results in light of this algorithm.

More precisely, we have generalised an integrated EM algorithm by preserving its monotonicity property for a wide range of hyperparameters. A particularly striking fact is that the monotonicity property holds for $\alpha \in [0, 1)$, hence updates akin to an EM procedure can be derived past the traditional case of likelihood optimisation. As we shall see, the additional layers of flexibility obtained in Algorithm 1 will also have important practical consequences due to the underlying gradient-based structure behind the mixture weights (and the mixture components parameters updates in the GMMs case) we have uncovered.

Table 1 summarises the main improvements of our framework compared to the existing literature and in the coming section, we revisit the maximisation and the gradient-based approaches when the variational family is based on the exponential family.

# 6. Exponential Family Distributions: a Closer Look

In this section, we state generic results for the maximisation and the gradient-based approaches in the important case where the variational family is based on the exponential family. Those results will in particular enable us to show how the update formulas in Examples 3 and 4 are derived, before allowing us to incorporate novel variational families in our framework such as multivariate Student's $t$ densities. To do so, let us first review the exponential family and some of its well-known properties.

## 6.1 Notation and Useful Definitions

We start with the definition of an exponential family distribution.
<table><caption>Table 1: Key improvements of our framework for mixture models optimisation compared to related works (prev. stands for previously in the literature).</caption><thead><tr><th></th><th>Improvements of our framework</th></tr></thead><tbody><tr><td>Rényi Gradient Descent<br>Li and Turner (2016)</td><td>Simultaneous optimisation w.r.t. (λ, Θ)<br>(prev. mixture weights λ optimisation not considered)<br>For GMMs : maximisation approach encompasses Rényi<br>Gradient Descent and provides covariance matrices updates</td></tr><tr><td>Power Descent<br>Daudel et al. (2021)</td><td>Simultaneous optimisation w.r.t. (λ, Θ)<br>(prev. (Θ<sub>n</sub>)<sub>n≥1</sub> constant)</td></tr><tr><td>M-PMC algorithm<br>Cappé et al. (2008)</td><td>Extension of an Integrated EM algorithm to:<br>α ∈ [0, 1], η<sub>n</sub> ∈ (0, 1], (α - 1)κ<sub>n</sub> ≥ 0 and b<sub>j,n</sub> ≥ 0<br>(prev. α = 0, η<sub>n</sub> = 1, κ<sub>n</sub> = 0 and b<sub>j,n</sub> = 0)</td></tr></tbody></table>

**Definition 1 (Exponential family)** Let $\nu$ be a $\sigma$-finite measure on $(\mathcal{Y}, \mathcal{Y})$, $E$ be a Euclidean space endowed with an inner product $\langle \cdot, \cdot \rangle_E$ and its corresponding norm $\|\cdot\|_E$, $h$ be a non-negative function defined on $\mathcal{Y}$ and $S: \mathcal{Y} \to E$ be a measurable function. In its canonical form, a member of the exponential family is defined by the parametric probability density

$$k^{(o)}(\zeta, y) = h(y) \exp (\langle \zeta, S(y) \rangle_E - A(\zeta)), \quad y \in \mathcal{Y}, \zeta \in E_0$$

where

$$E_0 := \left\{ \zeta \in E : \int_{\mathcal{Y}} h(y) e^{\langle \zeta, S(y) \rangle_E} \nu(dy) < \infty \right\}$$

and where the values taken by $A$ on the subset $E_0$ are obtained from the normalising constraint $\int k^{(o)}(\zeta, \cdot) d\nu = 1$. In this setting, $\nu$ is called the dominating measure and $S$ the natural statistic. In what follows, $\text{Int}E_0$ denotes the interior of $E_0$.

Recall now that $S$ defines an exhaustive statistic for this model, that $E_0$ is a convex subset of $E$ and that $A$ is infinitely differentiable in $\text{Int}E_0$. Moreover, for all $\zeta \in \text{Int}E_0$, the expectation and covariance operator of $S$ under the density $k^{(o)}(\zeta, \cdot)$ are the gradient $\nabla A$ and the Hessian $\nabla^T A$ of $A$ taken at $\zeta$, respectively. For later reference, we also recall

$$\nabla A(\zeta) = \int_{\mathcal{Y}} S(y) k^{(o)}(\zeta, y) \nu(dy) \quad (32)$$

and we introduce the following assumption.

(B1) The kernel density $k^{(o)}$ defined on $E_0 \times \mathcal{Y}$ is a member of the exponential family as in Definition 1. Moreover, we assume that:

(a) The function $h$ on $\mathcal{Y}$ is positive.

(b) There is no affine hyperplane of $E$ to which $S(y)$ belongs for $\nu$-almost all $y \in \mathcal{Y}$.

Notice that the assumption (B1)(a) can be circumvented by changing the dominating measure $\nu$ and therefore the assumptions (B1)(a) and (B1)(b) should be considered together. These
assumptions notably imply that the covariance operator of S under the density $k^{(o)}(\zeta, \cdot)$ is positive definite and thus so is $\nabla\nabla^T A(\zeta)$ for all $\zeta \in \text{Int} E_0$. In addition, one often relies on a non-canonical parameterisation (see for example Example 1), leading to the kernel density

$$k(\theta, y) = k^{(o)}(v(\theta), y) = h(y) \exp \left( \langle v(\theta), S(y) \rangle_E - A \circ v(\theta) \right), \quad \theta \in \mathcal{T}, y \in \mathcal{Y} \quad (33)$$

where $v$ maps the parameter space $\mathcal{T}$ to a subset of $\text{Int} E_0$. In that case, the assumption (B1)(a) will ensure that the assumption made on $k$ in (A1) holds.

## 6.2 The Maximisation Approach for the Exponential Family

Recall that the maximisation approach proposed in Corollary 1 aims at solving the optimisation problem (7) given by

$$\theta_{n+1} = \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_{\mathcal{Y}} \left[ \varphi_n^{(\alpha)}(y) + b_n k(\theta_n, y) \right] \log \left( \frac{k(\theta, y)}{k(\theta_n, y)} \right) \nu(dy), \quad n \ge 1.$$

Furthermore, in the more general case of mixture models described in Corollary 4, the maximisation approach leads to the argmax problem (25). This second argmax problem is similar to (7) in the sense that it updates $\Theta_n = (\theta_{j,n})_{1 \le j \le J}$ by solving J component-wise argmax problems of the same form as (7) (with $\theta_{n+1}, \varphi_n^{(\alpha)}, b_n$ and $\theta_n$ being replaced by $\theta_{j,n+1}, \varphi_{j,n}^{(\alpha)}, b_{j,n}$ and $\theta_{j,n}$, for $j = 1...J$):

$$\theta_{j,n+1} = \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_{\mathcal{Y}} \left[ \varphi_{j,n}^{(\alpha)}(y) + b_{j,n} k(\theta_{j,n}, y) \right] \log \left( \frac{k(\theta, y)}{k(\theta_{j,n}, y)} \right) \nu(dy), \quad j = 1...J, n \ge 1.$$

Our goal is now to solve the argmax problems above for a member of the exponential family. To do so, we first state a theorem in which the kernel $k$ is in its canonical form, that is $v$ is the identity mapping in (33) and $k(\theta, \cdot) = k^{(o)}(\zeta, \cdot)$. Other parameterisations will then be built using this theorem and as a result we will deduce corollaries that include non-canonical parameterisations and that are applicable to the general case of mixture models. All the proofs of the results stated in Section 6.2 are deferred to Appendix C.2.

**Theorem 4** Let $k^{(o)}$ satisfy (B1) and $\mathcal{O}$ be an open subset of $E_0$. Let $\check{\varphi}$ be a probability density function with respect to the measure $\nu$ such that

$$\int \|S\|_E \check{\varphi} \, d\nu < \infty. \quad (34)$$

Let $\zeta_0 \in \mathcal{O}$ and $b \ge 0$. Then the argmax problem

$$\underset{\zeta \in \mathcal{O}}{\operatorname{argmax}} \int_{\mathcal{Y}} \left[ \check{\varphi}(y) + b k^{(o)}(\zeta_0, y) \right] \log \left( \frac{k^{(o)}(\zeta, y)}{k^{(o)}(\zeta_0, y)} \right) \nu(dy) \quad (35)$$

admits at least a solution if and only if

$$\mathbf{s}^* := \frac{1}{1+b} \int S \check{\varphi} \, d\nu + \frac{b}{1+b} \nabla A(\zeta_0) \quad (36)$$

belongs to the image set $\nabla A(\mathcal{O})$, in which case (35) has a unique solution $\zeta^*$ defined by

$$\nabla A(\zeta^*) = \mathbf{s}^*. \quad (37)$$
Theorem 4 shows the equivalence between the argmax problem (35) and the equation (37) under the canonical parameterisation. From there, using the expression of $\nabla A$ given in (32) and considering the non-canonical exponential family probability density (33), we can interpret the equation (37) as an equality between two means of the statistic $S$ computed under two different distributions. Setting $\zeta^* = v(\theta^*)$ and $\zeta_0 = v(\theta_0)$, we indeed obtain

$$ \int_Y S(y)k(\theta^*, y)v(dy) = \int_Y S(y)\psi(y)v(dy) \quad (38) $$

where $\psi$ is the mixture

$$ \psi(y) = \frac{1}{1+b}\check{\varphi}(y) + \frac{b}{1+b}k(\theta_0, y), \quad y \in \mathcal{Y}. $$

As exemplified in the following corollary, an adequate parameterisation $\zeta = v(\theta)$ may then lead to a simple solution of (38), which in turn provides a solution to the argmax problem

$$ \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_{\mathcal{Y}} [\check{\varphi}(y) + bk(\theta_0, y)] \log \left( \frac{k(\theta, y)}{k(\theta_0, y)} \right) v(dy). \quad (39) $$

**Corollary 6 (Gaussian density)** Let $d \ge 1$ and $\mathcal{Y} = \mathbb{R}^d$. Let $\mathcal{M}_{>0}(d)$ denote the set of symmetric positive definite $d \times d$ matrices. We consider the case of a $d$-dimensional Gaussian density with $k(\theta, y) = \mathcal{N}(y; m, \Sigma)$, where $\theta = (m, \Sigma) \in \mathcal{T} = \mathbb{R}^d \times \mathcal{M}_{>0}(d)$. Let $\check{\varphi}$ be a probability density function on $\mathcal{Y}$ such that (34) holds, that is in the Gaussian case

$$ \int_{\mathcal{Y}} \|y\|^2 \check{\varphi}(y) \, dy < \infty. $$

Let $\theta_0 = (m_0, \Sigma_0) \in \mathcal{T}$ and $b \ge 0$. If $b = 0$, assume moreover that $\check{\varphi}$ is non-degenerate in the sense that its covariance matrix is positive definite. Then (39) has a unique solution $\theta^* = (m^*, \Sigma^*) \in \mathcal{T}$ defined by $m^* = \mathbb{E}[Y]$ and $\Sigma^* = \operatorname{Cov}(Y)$, where $Y$ is a random variable valued in $\mathbb{R}^d$ with density

$$ \psi(y) = \frac{1}{1+b}\check{\varphi}(y) + \frac{b}{1+b}\mathcal{N}(y; m_0, \Sigma_0), \quad y \in \mathcal{Y}. $$

It now remains to choose $\check{\varphi}$ adequately to relate (39) to (7), or rather to its generalisation (25). That way, we will in particular get that the update formulas given in Examples 1 and 3 are straightforward consequences of Corollary 6. This is the purpose of Corollary 7 below, in which we also rewrite the assumption (34) made on $\check{\varphi}$ under a more convenient form.

**Corollary 7** Let $k^{(0)}$ satisfy (B1). Further assume that the kernel $k$ is of the form (33) with a one-to-one mapping $v$ defined on $\mathcal{T}$ such that its image $v(\mathcal{T})$ is an open subset of $E_0$. Let $\alpha \in [0, 1)$ and let $p: \mathcal{Y} \to \mathbb{R}_+$ be such that

$$ 0 < \int_{\mathcal{Y}} \left( 1 + \|S(y)\|_{E}^{1/(1-\alpha)} \right) p(y) v(dy) < \infty \quad (40) $$

holds. Let $J \in \mathbb{N}$, $n \ge 1$ and $(\lambda_n, \Theta_n) \in S_J^+ \times \mathcal{T}^J$. Then, for all $j = 1 \dots J$, (34) holds with $\check{\varphi} = \check{\varphi}_{j,n}^{(\alpha)}$ as defined in (16) and setting $\gamma_{j,n} = \int \varphi_{j,n}^{(\alpha)} d\nu / (\int \varphi_{j,n}^{(\alpha)} d\nu + b_{j,n})$, we have that

$$ \underset{\theta \in \mathcal{T}}{\operatorname{argmax}} \int_{\mathcal{Y}} \left[ \varphi_{j,n}^{(\alpha)}(y) + b_{j,n}k(\theta_{j,n}, y) \right] \log \left( \frac{k(\theta, y)}{k(\theta_{j,n}, y)} \right) v(dy) \quad (41) $$
admits at least a solution if and only if

$$ \mathbf{s}_j^* := \gamma_{j,n} \int S \tilde{\varphi}_{j,n}^{(\alpha)} \, d\nu + (1 - \gamma_{j,n}) \nabla A \circ v(\theta_{j,n}) \quad (42) $$

belongs to the image set $\nabla A \circ v(\mathbf{T})$, in which case (41) has a unique solution $\theta_j^*$ defined by

$$ \nabla A \circ v(\theta_j^*) = \mathbf{s}_j^* $$

and we can set $\theta_{j,n+1} = \theta_j^*$ for all $j = 1 \dots J$ in (25) of Corollary 4.

Since the argmax problem (41) is immune to change of variables thanks to its argmax form, the maximisation approach can be solved using the canonical parameter $\zeta$ in Theorem 4 and transported to the parameter $\theta$ via the one-to-one mapping $v$ as written in Corollary 7.

The updates from Examples 1 and 3 then follow from Corollary 7, paired up with Corollary 6 and with the fact that (40) simplifies to (8) (since in the Gaussian case we have set $S(y) = (y, -yy^T/2)$ for all $y \in \mathbf{Y}$). More generally, we can always solve the argmax problem (41), that is the argmax problem of Corollary 4 with $k$ belonging to the exponential family, for $b_{j,n} > 0$ large enough.<fnref n="1" />

We now move on to the study of the gradient-based approach for the exponential family.

## 6.3 The Gradient-based Approach for an Exponential Family Distribution

Remember that to construct the sequence $(\theta_n)_{n \ge 1}$ according to the gradient-based approach from Corollary 2 we need to compute the gradient of the function $g_n$ defined in (11) by

$$ g_n(\theta) = \int_{\mathbf{Y}} \frac{\varphi_n^{(\alpha)}(y)}{\alpha - 1} \log \left( \frac{k(\theta, y)}{k(\theta_n, y)} \right) \nu(dy), \quad \theta \in \mathbf{T}, $$

and verify that $g_n$ satisfies a $\beta_n$-smoothness condition, which leads to updates of the form

$$ \theta_{n+1} = \theta_n - \frac{\gamma_n}{\beta_n} \nabla g_n(\theta_n), \quad n \ge 1 $$

with $\gamma_n \in (0, 1]$. Furthermore, the extension to mixture models in Corollary 5 involves the gradient of $J$ functions $g_{j,n}$ that resemble $g_n$ ($g_n, \varphi_n^{(\alpha)}$ and $\theta_n$ are replaced by $g_{j,n}, \varphi_{j,n}^{(\alpha)}$ and $\theta_{j,n}$ in the definition of $g_n$ above) and where each function $g_{j,n}$ is assumed to be $\beta_{j,n}$-smooth.

Let us now solve the gradient-based approach for a member of the exponential family. Like in the previous section, we start with a result handling the case where the kernel $k$ is in its canonical form (that is, $v$ is the identity mapping in (33) and $k(\theta, \cdot) = k^{(o)}(\zeta, \cdot)$). All the proofs of the results stated in Section 6.3 are deferred to Appendix C.3.

**Proposition 3** Let $k^{(o)}$ satisfy (B1). Let $\tilde{\varphi}$ be a probability density function with respect to $v$ satisfying (34). Let $\zeta_0 \in \text{Int} E_0$ and set

$$ g^{(o)}(\zeta) = - \int_{\mathbf{Y}} \tilde{\varphi}(y) \log \left( \frac{k^{(o)}(\zeta, y)}{k^{(o)}(\zeta_0, y)} \right) \nu(dy), \quad \zeta \in \text{Int} E_0. \quad (43) $$

<footnote n="1">1. Under the assumptions of Corollary 7, $v(\mathbf{T})$ is an open subset of $E_0$ and $\nabla A$ is a $C^\infty$-diffeomorphism on $\text{Int} E_0$. Hence, $\nabla A \circ v(\mathbf{T})$ is an open subset of $E$. Since $\mathbf{s}_j^*$ in (42) tends to $\nabla A \circ v(\theta_{j,n})$ as $b_{j,n} \to \infty$, $\mathbf{s}_j^*$ ends up belonging to $\nabla A \circ v(\mathbf{T})$ for $b_{j,n}$ large enough. The desired result follows by using that the regularisation parameter $b_{j,n} > 0$ in Corollary 4 can be freely chosen at each step $n \ge 1$.</footnote>
Then, for any convex subset $C_0 \subseteq \text{Int}E_0$, the function $g^{(o)}$ is $\beta_0$-smooth over $C_0$ with

$$ \nabla g^{(o)}(\zeta) = \nabla A(\zeta) - \int S \varphi \, d\nu \quad \text{and} \quad \beta_0 \ge \sup_{\zeta \in C_0} \|\nabla \nabla^T A(\zeta)\|_{\text{op}} \quad (44) $$

where $\|\cdot\|_{\text{op}}$ denotes the operator norm.

Letting $C_0$ be the line segment $\{\zeta_0 - \frac{t}{\beta_0}(\nabla A(\zeta_0) - \int S \varphi \, d\nu) : t \in [0, 1]\}$ in Proposition 3, we have that $C_0$ stays in $\text{Int}E_0$ for $\beta_0$ large enough. Since the Hessian of $A$ is continuous on $\text{Int}E_0$, the inequality in (44) can also be satisfied. Hence, the following gradient step

$$ \zeta = \zeta_0 - \frac{\gamma}{\beta_0} \left( \nabla A(\zeta_0) - \int_Y S \varphi \, d\nu \right) \quad (45) $$

can always be performed for all $\gamma \in (0, 1]$ and we can always apply the gradient-based approach when the kernel $k$ is in its canonical form. When considering non-canonical exponential family probability densities of the form (33) and under some reasonable assumptions on the mapping $v$ (for instance, if $v$ is a $C^2$ diffeomorphism), one is left with choosing between two alternative ways of exploiting Proposition 3:

(i) Apply a gradient step on the canonical parameter $\zeta$ using the objective $g^{(o)}$ in (43) and translate it into an update of the parameter $\theta$ via a change of variable. In this case, we set $\zeta_0 = v(\theta_0)$, we apply the step (45) and we then perform the change of variable $\theta = v^{-1}(\zeta)$ to get back to a parameter in $\mathbb{T}$. Overall, this leads to an update of the form

$$ \theta = v^{-1} \left( v(\theta_0) - \frac{\gamma}{\beta_0} \left( \nabla A \circ v(\theta_0) - \int_Y S \varphi \, d\nu \right) \right). $$

(ii) Apply a gradient step directly on the parameter $\theta$ using the objective

$$ g(\theta) = - \int \varphi(y) \log \left( \frac{k(\theta, y)}{k(\theta_0, y)} \right) dy, \quad \theta \in \mathbb{T}. $$

In that case, since $g = g^{(o)} \circ v$ where $g^{(o)}$ is again as in (43) with $\zeta_0 = v(\theta_0)$, this leads to an update of the form

$$ \theta = \theta_0 - \frac{\gamma}{\beta} \nabla v(\theta_0) \cdot \left( \nabla A \circ v(\theta_0) - \int_Y S \varphi \, d\nu \right). $$

In both approaches, $\gamma \in (0, 1]$ and the smoothness indices $\beta_0$ and $\beta$ have to be taken large enough in order to guarantee a decrease of the objective function. In practice, it is not clear which parameterisation leads to the simpler or more efficient algorithm. We investigate a particular setting in the following result.

**Corollary 8 (Gaussian density with known covariance matrix)** Let $d \ge 1$ and $Y = \mathbb{R}^d$. We consider the exponential family of a $d$-dimensional Gaussian density with $k(\theta, y) = N(y; \theta, \Sigma)$, where $\theta \in \mathbb{T} := \mathbb{R}^d$ and $\Sigma$ is a (known) covariance matrix in $M_{>0}(d)$. Let $\varphi$ be a probability density function w.r.t. $v$ such that (34) holds, that is in the Gaussian case

$$ \int \|y\| \varphi(y) dy < \infty. $$
Set $v(\theta) = \Sigma^{-1}\theta$ and define the canonical kernel $k^{(o)}$ with $h(y) = (2\pi)^{-d/2} |\Sigma|^{-1/2} e^{-y^T\Sigma^{-1}y/2}$, $S(y) = y$ and $A(\zeta) = \zeta^T\Sigma\zeta/2$. Then (33) holds and the methods (i) and (ii) lead to the following gradient-based updates, respectively,

$$ \theta = \theta_0 - \frac{\gamma}{\beta_0} \Sigma \left( \theta_0 - \int y \tilde{\varphi}(y) \, dy \right) \qquad (46) $$

$$ \theta = \theta_0 - \frac{\gamma}{\beta} \Sigma^{-1} \left( \theta_0 - \int y \tilde{\varphi}(y) \, dy \right) \qquad (47) $$

where $\beta_0$ is the largest eigenvalue of $\Sigma$ and $\beta$ is the largest eigenvalue of $\Sigma^{-1}$. They correspond to the smoothness indices of $g^{(o)}$ and $g$ over $\mathbb{R}^d$, respectively.

Corollary 8 illustrates how, even in a simple framework, the gradient-based updates strongly depend on the parameterisation chosen. They in fact coincide only if $\Sigma$ is scalar in the above corollary, in which case $\Sigma/\beta_0$ and $\Sigma^{-1}/\beta$ above are both equal to the identity matrix.

Following the reasoning of Section 6.2, we can then apply Corollary 8 with $\tilde{\varphi} = \tilde{\varphi}_n^{(\alpha)}$ (and more generally with $\tilde{\varphi} = \tilde{\varphi}_{j,n}^{(\alpha)}$). As (34) is implied by (13) (using Lemma 8 in Appendix C.2 and that here $S(y) = y$), this enables us to deduce the updates in Examples 2 and 4. Note that we did not introduce a convex subset $C_0$ in Corollary 8 (this is due to the fact that the smoothness index is constant over the whole parameter space for both gradient-based updates in that case), which is why $C_0$ does not appear in Examples 2 and 4.

Let us next put into perspective Corollary 6, in which we considered a Gaussian density with varying mean vector and covariance matrix, and Corollary 8, where the covariance matrix is assumed to be known. While Corollary 8 kept the computations straightforward and provided gradient-based updates that do not require to introduce a convex subset $C_0$, this is no longer the case if we take the same exact setting as in Corollary 6. Indeed, when $\theta = (m, \Sigma) \in \mathbb{T} = \mathbb{R}^d \times \mathcal{M}_{>0}(d)$ and under the condition (34), the gradient of $g$ is given by

$$ \nabla g(\theta) = \left( \Sigma^{-1} \left( m - \int y \tilde{\varphi}(y) \, dy \right), \frac{1}{2} \Sigma^{-1} \left( \int y y^T \tilde{\varphi}(y) \, dy - \Sigma + mm^T \right) \Sigma^{-1} \right). $$

The smoothness index $\beta$ is now no longer constant here, which makes it much more involved to choose a convenient convex subset $C_0$ and to bound the smoothness index over it.

At this stage, we can solve the maximisation and gradient-based approaches when the kernel $k$ belongs to the exponential family for variational families as large as (finite) mixture models. The maximisation approach appears to be preferable to the gradient-based one as it does not require a smoothness assumption and does not depend on the parameterisation.

In the next section, we investigate how our maximisation approach may further generalise to include extensions of the exponential family distribution such as Student's $t$ distributions, hence further motivating the maximisation approach over the gradient-based approach.

## 6.4 Extension to Linear Mixture (LM) Models

The goal of this section is to examine what becomes of the maximisation approach for a specific extension of the exponential family distribution, which we will refer to as exponential linear mixture (ELM) family. The ELM family will in particular encompass the Student's
t distribution with mean m, covariance matrix Σ and a degrees of freedom which can be defined as the continuous mixture

$$t(y; m, \Sigma, a) = \int_{0}^{\infty} \mathcal{N}(y; m, z^{-1}\Sigma) \tilde{\tau}_{a}(dz),$$

where $\tilde{\tau}_a$ denotes the $\chi^2$ distribution with a degrees of freedom. As a result, building a maximisation approach within the ELM family will lead to new iterative schemes for variational families that go beyond the cases considered thus far. To this end, let us first provide a precise definition of the ELM family we want to study and introduce the main assumptions of this section. All the proofs from Section 6.4 are deferred to Appendix C.4.

### 6.4.1 THE (E)LM FAMILY: DEFINITION AND NOTATION

Let $F$ be an Euclidean space endowed with the inner product $\langle \cdot, \cdot \rangle_F$ and norm $\| \cdot \|_F$. We denote by $\mathcal{L}(F, E)$ and $\mathcal{L}(F)$ the (finite-dimensional) linear spaces of linear operators from $F$ to $E$ and from $F$ to itself, respectively. We now provide the definition of the (E)LM family.

**Definition 2 ( (E)LM family) Let $k^{(o)}: E_0 \times \mathbb{Y} \to \mathbb{R}_+$ be a kernel density w.r.t. the dominating measure $\nu$, with $E_0 \subseteq E$ and let $\Upsilon: \tilde{\mathbb{T}} \to F$. Set $\mathcal{Z} = \mathcal{L}(F, E)$ and denote by $\mathcal{Z}$ its Borel $\sigma$-field. Let $M$ be a subset of probability measures on $(\mathcal{Z}, \mathcal{Z})$ satisfying**

(i) for all $\tau, \tau'$ in $M$, the distributions $\tau$ and $\tau'$ are equivalent;

(ii) for all $\tau \in M$ and $\tau$-almost all $\ell \in \mathcal{Z}$, the image set $\ell \circ \Upsilon(\tilde{\mathbb{T}})$ is included in $E_0$.

A member of the linear mixture (LM) family with canonical kernel $k^{(o)}$, mixing class $M$ and parameter mapping $\Upsilon$, is defined for a parameter $\theta = (\vartheta, \tau) \in \mathbb{T} = \tilde{\mathbb{T}} \times M$ by the density w.r.t. $\nu$

$$k^{(1)}(\theta, y) = \int_{\mathcal{Z}} k^{(o)}(\ell \circ \Upsilon(\vartheta), y) \tau(d\ell), \quad y \in \mathbb{Y}. \qquad (48)$$

If moreover $k^{(o)}$ satisfies (B1), we say that $k^{(1)}(\theta, \cdot)$ is a member of the exponential linear mixture (ELM) family.

Here, Assumption (i) implies that we only need to check Assumption (ii) for one $\tau$ in $M$. As for Assumption (ii), it ensures that for all $\vartheta \in \tilde{\mathbb{T}}$ and all $\tau \in M$, the function $(y, \ell) \mapsto k^{(o)}(\ell \circ \Upsilon(\vartheta), y)$ appearing in (48) is a probability density function w.r.t. $\nu \otimes \tau$, so that for all $\theta \in \mathbb{T}$, the function $y \mapsto k^{(1)}(\theta, y)$ is a probability density function w.r.t. $\nu$.

Furthermore, the mapping $\Upsilon$ in Definition 2 allows us to introduce various parameterisations. It is also important to note that for $\theta = (\vartheta, \tau) \in \mathbb{T}$, the probability density function $k^{(1)}(\theta, \cdot)$ only depends on $\xi = \Upsilon(\vartheta) \in \Upsilon(\tilde{\mathbb{T}})$ and $\tau \in M$. As we shall see later, Student's t distributions in particular will fit this general setting. Before that, let us state conditions leading to a systematic decrease in $\Psi_\alpha$ when the kernel $k^{(1)}$ is as in Definition 2 and from there, let us see how a maximisation approach for the ELM family can be derived.
### 6.4.2 MONOTONIC DECREASE CONDITIONS FOR THE (E)LM FAMILY

In the case of an ELM family, we are not able to directly solve the argmax problem (25) as we did for the exponential family in Section 6.2. Instead, given an (E)LM family $k^{(1)}$ as in Definition 2, we come back to a monotonic decrease condition of the form

$$ \int_Y \tilde{\varphi}(y) \log \left( \frac{k^{(1)}(\theta, y)}{k^{(1)}(\theta_0, y)} \right) \nu(dy) \ge 0. \quad (49) $$

As detailed in the remark below, for selected choices of $\tilde{\varphi}$ and of $(\theta, \theta_0)$, (49) can indeed be linked to the conditions (5) and (19) appearing in Theorem 1 and Theorem 2 respectively.

**Remark 1** Let $\alpha \in [0, 1)$. Taking $\tilde{\varphi} = \tilde{\varphi}_n^{(\alpha)}$ defined by (6) and $(\theta, \theta_0) = (\theta_{n+1}, \theta_n)$, (49) becomes the condition (5) of Theorem 1. Applying (49) for $j = 1, \dots, J$ with $\tilde{\varphi} = \tilde{\varphi}_{j,n}^{(\alpha)}$ defined by (16) and $(\theta, \theta_0) = (\theta_{j,n+1}, \theta_{j,n})$ yields the condition (19) of Theorem 2.

Let $D_{\tau', \tau} = \frac{d\tau'}{d\tau}$ denote the Radon-Nikodym derivative of $\tau'$ w.r.t. $\tau$. We then have the following proposition (we defer the proofs of the results from Section 6.4.2 to Appendix C.4.1).

**Proposition 4 (Conditions for a monotonic decrease within the (E)LM family)**

Let $k^{(1)}$ be an LM family as in Definition 2. Let $\tilde{\varphi}: Y \to \mathbb{R}_+$ be a probability density function with respect to $\nu$. Let $\vartheta_0, \vartheta \in \tilde{T}$ and $\tau_0, \tau \in M$ such that

$$ \int_Y \left( \int_Z \varphi(\ell, y) \log (D_{\tau, \tau_0}(\ell)) \tau_0(d\ell) \right) \nu(dy) \ge 0, \quad (50) $$

$$ \int_Y \left( \int_Z \varphi(\ell, y) \log \left( \frac{k^{(0)}(\ell \circ \Upsilon(\vartheta), y)}{k^{(0)}(\ell \circ \Upsilon(\vartheta_0), y)} \right) \tau_0(d\ell) \right) \nu(dy) \ge 0 \quad (51) $$

where we set

$$ \varphi(\ell, y) = \tilde{\varphi}(y) \frac{k^{(0)}(\ell \circ \Upsilon(\vartheta_0), y)}{k^{(1)}(\theta_0, y)}. \quad (52) $$

Then (49) holds with $\theta = (\vartheta, \tau)$ and $\theta_0 = (\vartheta_0, \tau_0)$.

One possible way to obtain (50) is to observe that since the left-hand side of (50) is zero for $\tau = \tau_0$, (50) is fulfilled by setting

$$ \tau = \underset{\tau' \in M}{\operatorname{argmax}} \int_Y \left( \int_Z \varphi(\ell, y) \log (D_{\tau', \tau_0}(\ell)) \tau_0(d\ell) \right) \nu(dy), \quad (53) $$

assuming that this argmax is well-defined. This will notably be done later in order to get the update formula (70) of Example 5 (see the proof of Corollary 10 in Appendix C.4).

As for the updating of $\vartheta$, we can again adopt a maximisation approach to derive a convenient $\vartheta$. This is the purpose of the following result, where we directly update $\xi = \Upsilon(\vartheta)$ since $\Upsilon$ is a known mapping in Definition 2.
**Corollary 9** Consider an LM family as in Definition 2 and $\theta_0 = (\vartheta_0, \tau_0) \in \tilde{T} \times M$. Set $\xi_0 = \Upsilon(\vartheta_0)$, let $\tilde{\varphi}: Y \to \mathbb{R}_+$ be a probability density function with respect to $\nu$ and define $\varphi: Z \times Y \to \mathbb{R}_+$ by (52). Let $\vartheta \in \tilde{T}$ be such that $\Upsilon(\vartheta)$ is a solution of the argmax problem

$$ \underset{\xi \in \Upsilon(\tilde{T})}{\operatorname{argmax}} \int_Y \left( \int_Z \left[ \varphi(\ell, y) + b_0(\ell) k^{(o)}(\ell(\xi_0), y) \right] \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \tau_0(d\ell) \right) \nu(dy), \quad (54) $$

where $b_0$ is any non-negative function defined on $Z$ such that $\int b_0 d\tau_0 < \infty$. Then (51) holds.

As we will see, (53) and (54) can be achieved for continuous mixtures such as the one used to define Student's $t$ distributions but they require a thorough analysis. To this end, we next formulate and solve a maximisation approach for a kernel $k^{(1)}$ belonging to the ELM family.

### 6.4.3 THE MAXIMISATION APPROACH FOR THE ELM FAMILY

Let us show that we can solve an argmax of the form (54) for an adequate choice of $b_0$ by fully exploiting the fact that $k^{(o)}(\zeta, y)$ is a canonical kernel of the exponential family distribution which satisfies (B1) and by relying on additional assumptions to be introduced alongside some helpful notation. Our first assumption is the following.

(C1) The image set $\Upsilon(\tilde{T})$ is an open subset of $F$ and for all $\tau \in M$, $\tau$-a.e. $\ell$ is full rank.

(C1) guarantees that, for all $\vartheta \in \tilde{T}$, $\tau \in M$ and $\tau$-almost all $\ell$ in $Z$, $\ell \circ \Upsilon(\vartheta)$ belongs to $\mathrm{Int}E_0$, the interior set of $E_0$ (see Lemma 9 in Appendix C.4 for details). This will come in handy in the upcoming derivations and we now introduce some helpful notation, before presenting our next two assumptions. We let $\|\cdot\|_{op}$ denote operator norms, for instance, for any $\ell \in Z$,

$$ \|\ell\|_{op} = \sup \{ \|\ell(x)\|_F : x \in E, \|x\|_E \le 1 \}. $$

We further denote the adjoint of the linear operator $\ell$ by $\ell^T$, for instance, in the case $\ell \in Z$, $\ell^T \in \mathcal{L}(E, F)$ is defined by

$$ \langle \ell(x), y \rangle_E = \langle x, \ell^T(y) \rangle_F \quad \text{for all } (x, y) \in F \times E \quad (55) $$

and, for convenience, we also denote

$$ \ell^{\nabla} := \ell^T \circ \nabla A \circ \ell. \quad (56) $$

Here, $A$ is the function appearing in the definition of an exponential family (Definition 1). Recall that $A$ is infinitely differentiable on $\mathrm{Int}E_0$ so that by (C1), $\ell^{\nabla}$ is well-defined on $\Upsilon(\tilde{T})$ for $\tau$-a.e. $\ell$ with $\tau \in M$. Furthermore, we use $\mathbb{P}_{\vartheta, \tau}$ and $\mathbb{E}_{\vartheta, \tau}$ to denote the probability and its corresponding expectation under which the pair $(Y, L)$ has density $(y, \ell) \mapsto k^{(o)}(\ell \circ \Upsilon(\vartheta), y)$ with respect to $\nu \otimes \tau$. Namely, for any non-negative measurable $g$ on $Y \times Z$,

$$ \mathbb{E}_{\vartheta, \tau}[g(Y, L)] = \int g(y, \ell) k^{(o)}(\ell \circ \Upsilon(\vartheta), y) \nu(dy) \tau(d\ell). \quad (57) $$

We now introduce, for any $(\vartheta_0, \tau_0) \in \tilde{T} \times M$, two auxiliary functions $\tilde{m}_{\vartheta_0, \tau_0}$ and $\tilde{s}_{\vartheta_0, \tau_0}$, both defined from $Y$ to $\mathbb{R}_+$, such that the two following conditions hold.
(C2) For all $\vartheta_0 \in \tilde{T}$ and $\tau_0 \in M$, we have

$$ \mathbb{E}_{\vartheta_0, \tau_0} \left[ \|L\|_{\text{op}} \mid Y \right] \le \tilde{m}_{\vartheta_0, \tau_0}(Y) \quad \mathbb{P}_{\vartheta_0, \tau_0} - \text{a.s.} $$

(C3) For all $\vartheta_0 \in \tilde{T}$ and $\tau_0 \in M$, we have

$$ \forall \xi \in \Upsilon(\tilde{T}), \exists \epsilon, C > 0, \mathbb{E}_{\vartheta_0, \tau_0} \left[ e^{\epsilon \|L^{\nabla}(\xi)\|_F} \mid Y \right] \le C e^{\tilde{s}_{\vartheta_0, \tau_0}(Y)} \quad \mathbb{P}_{\vartheta_0, \tau_0} - \text{a.s.} $$

(C2) and (C3) are made to ensure that the integrals used to solve the argmax problem from Corollary 9 will be well-defined. We finally present our last assumption, which corresponds to some sort of identifiability assumption (see Lemma 10 in Appendix C.4 for details). It is used to obtain the uniqueness for the argmax problem we will solve.

(C4) For all $\xi \ne \xi' \in \Upsilon(\tilde{T})$ and $\tau_0 \in M$, we have $\tau_0(\{\ell \in \mathbb{Z} : \ell(\xi) \ne \ell(\xi')\}) > 0$.

We can now state the following result, whose proof can be found in Appendix C.4.2.

**Theorem 5** Consider an ELM family as in Definition 2 and let $\theta_0 = (\vartheta_0, \tau_0) \in \tilde{T} \times M$. Assume (C1)-(C4) hold. Let $\tilde{\varphi}$ be a probability density function w.r.t. $\nu$ such that

$$ \int_{\Upsilon} \tilde{\varphi}(y) \left( \|S(y)\|_E \tilde{m}_{\vartheta_0, \tau_0}(y) + e^{\tilde{s}_{\vartheta_0, \tau_0}(y)} \right) \nu(dy) < \infty, \quad (58) $$

define $\varphi: \mathbb{Z} \times \mathbb{Y} \to \mathbb{R}_+$ by (52) and set for all $(\ell, y, \xi) \in \mathbb{Z} \times \mathbb{Y} \times \Upsilon(\tilde{T})$ such that $\ell \circ \Upsilon(\vartheta_0) \in \text{Int} E_0$

$$ \tilde{\varphi}(\ell) = \int_{\Upsilon} \varphi(\ell, y) \nu(dy), \quad (59) $$

$$ w_{\tau_0}^{\nabla}(\xi) = \int_{\mathbb{Z}} \tilde{\varphi}(\ell) \ell^{\nabla}(\xi) \tau_0(d\ell). \quad (60) $$

Then $\tilde{\varphi} < \infty$, $\tau_0$-a.s., and $w_{\tau_0}^{\nabla}$ is well-defined from $\Upsilon(\tilde{T})$ to $F$ and one-to-one on $\Upsilon(\tilde{T})$, hence bijective from $\Upsilon(\tilde{T})$ to its image $w_{\tau_0}^{\nabla} \circ \Upsilon(\tilde{T})$, and satisfies, for all $\vartheta \in \tilde{T}$,

$$ w_{\tau_0}^{\nabla} \circ \Upsilon(\vartheta) = \mathbb{E}_{\vartheta, \tau_0} [L^T \circ S(Y) \tilde{\varphi}(L)] . \quad (61) $$

Moreover, setting $\xi_0 = \Upsilon(\vartheta_0)$, for any $b \in \mathbb{R}_+$, the argmax problem

$$ \underset{\xi \in \Upsilon(\tilde{T})}{\operatorname{argmax}} \int_{\Upsilon} \left( \int_{\mathbb{Z}} \left[ \varphi(\ell, y) + b \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \right] \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \tau_0(d\ell) \right) \nu(dy), \quad (62) $$

has at least a solution if and only if

$$ \mathbf{w}^* = \frac{1}{1+b} \int \varphi(\ell, y) \ell^T \circ S(y) \nu(dy) \tau_0(d\ell) + \frac{b}{1+b} w_{\tau_0}^{\nabla}(\xi_0) \quad (63) $$

belongs to $w_{\tau_0}^{\nabla} \circ \Upsilon(\tilde{T})$, in which case (62) has a unique solution $\xi^*$ defined by

$$ w_{\tau_0}^{\nabla}(\xi^*) = \mathbf{w}^*. \quad (64) $$
By setting $\xi^* = \Upsilon(\vartheta^*)$ and $\xi_0 = \Upsilon(\vartheta_0)$ as well as using (61), (64) can be interpreted as

$$ \mathbb{E}_{\vartheta^*, \tau_0} [L^T \circ S(Y) \tilde{\varphi}(L)] = \mathbb{E}_{\vartheta_0, \tilde{\tau}_0} [L^T \circ S(\tilde{Y})], \quad (65) $$

where $\tilde{\tau}_0$ denotes the probability having density $\tilde{\varphi}$ with respect to $\tau_0$ and $\tilde{Y}$ denotes a random variable valued in $\mathcal{Y}$ such that, under $\mathbb{P}_{\vartheta_0, \tilde{\tau}_0}$, the pair $(\tilde{Y}, L)$ is distributed according to the mixture density

$$ (y, \ell) \mapsto \psi(y, \ell) = \frac{1}{1+b} \varphi(\ell, y) + \frac{b}{1+b} k^{(o)}(\ell \circ \Upsilon(\vartheta_0), y) \tilde{\varphi}(\ell) $$

with respect to $\nu \otimes \tau_0$. Interestingly, the second marginals of these two mixands are $\tilde{\tau}_0$. The distribution of $(\tilde{Y}, L)$ under $\mathbb{P}_{\vartheta_0, \tilde{\tau}_0}$ can thus equivalently be defined by saying that $L \sim \tilde{\tau}_0$ and the conditional distribution of $\tilde{Y}$ given $L = \ell$ is the mixture with density

$$ y \mapsto \psi[y | \ell] = \frac{1}{1+b} \frac{\varphi(\ell, y)}{\tilde{\varphi}(\ell)} + \frac{b}{1+b} k^{(o)}(\ell \circ \Upsilon(\vartheta_0), y). $$

Using the definition of $\tilde{\tau}_0$ on the left-hand side of (65) as well, the latter reads

$$ \mathbb{E}_{\vartheta^*, \tilde{\tau}_0} [L^T \circ S(Y)] = \mathbb{E}_{\vartheta_0, \tilde{\tau}_0} [L^T \circ S(\tilde{Y})]. \quad (66) $$

Hence, Theorem 5 along with (66) can be used to solve the argmax problem (54) in the case where (B1) and (C1)-(C4) hold and for $b_0(\ell) = b \tilde{\varphi}(\ell)$ with $b$ being a non-negative constant. From there, it only remains to find ways to check that the condition (58) on $\tilde{\varphi}$ holds. This can notably be done when $\tilde{\varphi}$ takes one of the two forms listed in Remark 1 (see Lemma 8 in Appendix C.2 for details) and we now provide an example of particular interest.

**Example 5 (Finite mixture of Student's t distributions)** Let $\nu$ be the Lebesgue measure on $\mathcal{Y} = \mathbb{R}^d$, with $d \ge 1$. Let $\alpha \in [0, 1)$ and $p: \mathbb{R}^d \to \mathbb{R}_+$ such that

$$ 0 < \int_{\mathbb{R}^d} (1 + \|y\|^2)^{1/(1-\alpha)} p(y) dy < \infty. \quad (67) $$

Let $J \ge 1$ and consider a J-mixture family $\mathcal{Q}$ as in (15), with $k(\theta, \cdot)$ being the density of the Student's t distribution with parameter $\theta = (m, \Sigma, a) \in \mathcal{T} = \mathbb{R}^d \times \mathcal{M}_{>0}(d) \times (0, \infty)$, given by

$$ k(\theta, y) = \frac{\Gamma((a+d)/2)}{\Gamma(a/2) (a\pi)^{d/2} |\Sigma|^{1/2}} \left(1 + \frac{1}{a}(y-m)^T \Sigma^{-1} (y-m)\right)^{-(a+d)/2}. $$

Define a sequence $(\mu_n k)_{n \ge 1}$ valued in $\mathcal{Q}$ where $\mu_n = \mu_{\lambda_n, \Theta_n}$ with $\lambda_n \in S_J^+$ and $\Theta_n = (m_{j,n}, \Sigma_{j,n}, a_{j,n})_{j=1...J} \in \mathcal{T}^J$ such that, for all $n \ge 1$,

(i) given $(\eta_n)_{n \ge 1}$ valued in $(0, 1]$, $(\kappa_n)_{n \ge 1}$ valued in $(-\infty, 0]$ and $\lambda_1 \in S_J^+$, $(\lambda_n)_{n \ge 2}$ is defined by the iterative formula (21);
(ii) given $(\gamma_{j,n})_{j=1,...,J, n \ge 1}$ valued in $(0,1]$ and $\Theta_1 \in \mathcal{T}^J$, $(\Theta_n)_{n \ge 2}$ is defined by

$$m_{j,n+1} = \left( \int_{\mathbb{R}^d \times \mathbb{R}_+} z \psi_{j,n}(z, y) \, dy \, \check{\tau}_{a_{j,n}}(dz) \right)^{-1} \int_{\mathbb{R}^d \times \mathbb{R}_+} y z \psi_{j,n}(z, y) \, dy \, \check{\tau}_{a_{j,n}}(dz), \quad (68)$$

$$\Sigma_{j,n+1} = \int_{\mathbb{R}^d \times \mathbb{R}_+} (y - m_{j,n+1})(y - m_{j,n+1})^T z \psi_{j,n}(z, y) \, dy \, \check{\tau}_{a_{j,n}}(dz), \quad (69)$$

$$a_{j,n+1} = 2 \kappa^{-1} \left( \int_{\mathbb{R}^d \times \mathbb{R}_+} (z - \ln(z)) \varphi_{j,n}(z, y) \, dy \, \check{\tau}_{a_{j,n}}(dz) \right), \quad (70)$$

where, defining $\check{\varphi}_{j,n}^{(\alpha)}$ as in (16) and denoting $\kappa^{-1} : \mathbb{R} \to (0, \infty)$ the inverse mapping of $\kappa(x) = \log(x) + \Gamma'(x)/\Gamma(x)$, we set

$$\varphi_{j,n}(z, y) = \check{\varphi}_{j,n}^{(\alpha)}(y) \frac{e^{-\frac{z}{2}(y-m_{j,n})^T \Sigma_{j,n}^{-1} (y-m_{j,n})}}{(2\pi/z)^{d/2} |\Sigma_{j,n}|^{1/2} k(\theta_{j,n}, y)}, \quad (71)$$

$$\psi_{j,n}(z, y) = \gamma_{j,n} \varphi_{j,n}(z, y) + (1 - \gamma_{j,n}) \frac{e^{-\frac{z}{2}(y-m_{j,n})^T \Sigma_{j,n}^{-1} (y-m_{j,n})}}{(2\pi/z)^{d/2} |\Sigma_{j,n}|^{1/2}} \left( \int \varphi_{j,n}(z, y) \, dy \right), \quad (72)$$

$$\check{\tau}_a(dz) = \mathbb{1}_{\mathbb{R}_+}(z) \frac{(a/2)^{a/2}}{\Gamma(a/2)} z^{a/2-1} e^{-za} dz \quad \text{for any } a > 0, \quad (73)$$

Here, $\check{\tau}_a$ in (73) is the $\chi^2$ distribution with a degrees of freedom and $\kappa^{-1}$ in (70) is well-defined (this follows from Lemma 13 in Appendix C.4.3).

We then have the following result, whose proof is postponed to Appendix C.4.3.

**Corollary 10** In the setting of Example 5, if $\Psi_\alpha(\mu_1 k) < \infty$, then, at any time $n \ge 1$, we have $\Psi_\alpha(\mu_{n+1} k) \le \Psi_\alpha(\mu_n k)$.

We have stated results that solve our maximisation and gradient-based approaches when the variational family is based on the exponential family. However, all the update formulas we have obtained so far involve intractable integrals. In the next section, we focus on the case of Gaussian Mixture Models optimisation seen in Examples 3 and 4 and we investigate how our framework can be used to derive tractable algorithms in that case.

## 7. Gaussian Mixture Models (GMMs) Optimisation

We consider $d$-dimensional Gaussian mixture densities with $k(\theta_j, y) = \mathcal{N}(y; m_j, \Sigma_j)$, where $\theta_j = (m_j, \Sigma_j) \in \mathcal{T}$ denotes the mean and covariance matrix of the $j$-th Gaussian component density. By Theorem 2, a valid update formula at time $n \ge 1$ for $(\lambda_n)_{n \ge 1}$ is

$$\lambda_{j,n+1} = \frac{\lambda_{j,n} \left[ \int_{\mathcal{Y}} \varphi_{j,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1)\kappa_n \right]^{\eta_n}}{\sum_{\ell=1}^{J} \lambda_{\ell,n} \left[ \int_{\mathcal{Y}} \varphi_{\ell,n}^{(\alpha)}(y) \nu(dy) + (\alpha - 1)\kappa_n \right]^{\eta_n}}, \quad j = 1 \dots J. \quad (74)$$

Furthermore, as underlined in Example 3, a valid update at time $n \ge 1$ for the means $(m_{j,n+1})_{1 \le j \le J}$ and the covariances matrices $(\Sigma_{j,n+1})_{1 \le j \le J}$ is given by

$$\forall j = 1 \dots J, \quad m_{j,n+1} = (1 - \gamma_{j,n})m_{j,n} + \gamma_{j,n}\hat{m}_{j,n} \quad (75)$$

$$\Sigma_{j,n+1} = (1 - \gamma_{j,n})\Sigma_{j,n} + \gamma_{j,n}\hat{\Sigma}_{j,n} + \gamma_{j,n}(1 - \gamma_{j,n}) (\hat{m}_{j,n} - m_{j,n}) (\hat{m}_{j,n} - m_{j,n})^T$$
where $\hat{m}_{j,n} = \int_Y y \tilde{\varphi}_{j,n}^{(\alpha)}(y) \nu(dy)$, $\hat{\Sigma}_{j,n} = \int_Y y y^T \tilde{\varphi}_{j,n}^{(\alpha)}(y) \nu(dy) - \hat{m}_{j,n} \hat{m}_{j,n}^T$ with $\tilde{\varphi}_{j,n}^{(\alpha)} = \varphi_{j,n}^{(\alpha)} / \int \varphi_{j,n}^{(\alpha)} \nu(dy)$ and $\gamma_{j,n} \in (0, 1]$. In addition, another possible update for the means is

$$m_{j,n+1} = m_{j,n} + \gamma_{j,n} \frac{\int_{\mathbb{Y}} \lambda_{j,n} \varphi_{j,n}^{(\alpha)}(y) (y - m_{j,n}) \nu(dy)}{\int_{\mathbb{Y}} \mu_n k(y)^{\alpha} p(y)^{1-\alpha} \nu(dy)}, \quad j = 1 \dots J. \quad (76)$$

By virtue of Example 4, this update can be used in lieu of the means update written in (75) with $\gamma_{j,n}$ still in $(0, 1]$. Since the update (76) is linked to Gradient Descent steps for Rényi's $\alpha$-divergence minimisation (see Example 4), we will refer to this update as the *Rényi Gradient Descent* (RGD) update in the following. One the other hand, the update on the means written in (75) will be referred to as the *Maximisation Gradient* (MG) update.

Notice that the updates (74), (75) and (76) all involve intractable integrals. However, we can resort to approximate update rules and take advantage of the fact that the integrals appearing in these update formulas are of the form

$$\int_{\mathbb{Y}} \varphi_{j,n}^{(\alpha)}(y) g(y) \nu(dy), \quad (77)$$

where $g$ is a well-chosen function defined on $\mathbb{Y}$ and where $\varphi_{j,n}^{(\alpha)}$ is defined in (16) by $\varphi_{j,n}^{(\alpha)}(y) = k(\theta_{j,n}, y) (\mu_n k(y)/p(y))^{\alpha-1}$ for all $y \in \mathbb{Y}$. Many choices are indeed possible to approximate (77) and for simplicity we restrict ourselves to typical Importance Sampling estimation. Looking at the definition of $\varphi_{j,n}^{(\alpha)}$, a first idea would be to sample $M$ i.i.d. random variables $(Y_{m,n}^{(j)})_{1 \le m \le M}$ according to $k(\theta_{j,n}, \cdot)$ for $j = 1 \dots J$ and to use the unbiased estimator of (77)

$$\frac{1}{M} \sum_{m=1}^{M} \left( \frac{\mu_n k(Y_{m,n}^{(j)})}{p(Y_{m,n}^{(j)})} \right)^{\alpha-1} g(Y_{m,n}^{(j)}).$$

As this sampler depends on $j$, this option becomes very computationally heavy as $J$ and $M$ increase due to the *sampling cost* ($J \times M$) and to the *evaluation cost* ($\mu_n k$ and $p$ need to be evaluated $M \times J$ times each, where $\mu_n k$ is a sum over $J$). To reduce this computational bottleneck, it is then preferable to consider sequences of samplers $(q_n)_{n \ge n}$ that do not depend on $j$ at time $n$. More specifically, we discuss two possibilities.

(i) **Best sampler at time $n$ (IS-n).** Approximated update rules can be obtained by sampling $M$ i.i.d. random variables $(Y_{m,n})_{1 \le m \le M}$ according to the best approximation of $p$ we have a time $n$, that is $q_n = \mu_n k$. Then, the samples are shared throughout the $J$ mixture weights and mixture components parameters updates. This eases the computational burden, with both a sampling and an evaluation costs of $M$, as an unbiased estimator of (77) is

$$\frac{1}{M} \sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n}) g(Y_{m,n}) \quad \text{where} \quad \hat{\varphi}_{j,n}^{(\alpha)}(y) = \frac{k(\theta_{j,n}, y)}{q_n(y)} \left( \frac{\mu_n k(y)}{p(y)} \right)^{\alpha-1}, \quad y \in \mathbb{Y}. \quad (78)$$

(ii) **Uniform sampler (IS-unif).** Approximated update rules can be obtained by sampling according to the uniform sampler $q_n = J^{-1} \sum_{j=1}^{J} k(\theta_{j,n}, \cdot)$ and using the unbiased estimator (78). That way, the sampling and evaluation costs are also $M$ each. Contrary to IS-n, this sampler ensures a fair sampling among all components, as it entails sampling an index vector $[i_1, i_2, \dots, i_M]$ uniformly from $\{1, \dots, J\}$ whereas IS-n samples an index vector $[i_1, i_2, \dots, i_M]$ according to the mixture weights $\lambda_n$.
**Algorithm 3:** GMMs optimisation with the IS-n/IS-unif sampler

At iteration n,

1. Draw independently M samples $(Y_{m,n})_{1\le m\le M}$ from the proposal $q_n$.

2. For all $j = 1...J$,

(a) Set $\hat{m}_{j,n} = \frac{\sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n}) \cdot Y_{m,n}}{\sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n})}$ and $\hat{\Sigma}_{j,n} = \frac{\sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n}) \cdot Y_{m,n} Y_{m,n}^{T}}{\sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n})} - \hat{m}_{j,n} \hat{m}_{j,n}^{T}$.

(b) Choose one between the (MG) and (RGD) approaches

$$ (MG) \quad m_{j,n+1} = (1 - \gamma_n) m_{j,n} + \gamma_n \hat{m}_{j,n} $$

$$ (RGD) \quad m_{j,n+1} = m_{j,n} + \gamma_n \frac{\lambda_{j,n} \sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n}) \cdot (Y_{m,n} - m_{j,n})}{\sum_{j=1}^{J} \sum_{m=1}^{M} \lambda_{j,n} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n})} $$

and set

$$ \lambda_{j,n+1} = \frac{\lambda_{j,n} \left[ \sum_{m=1}^{M} \hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n}) + (\alpha - 1)\kappa_n \right]^{\eta_n}}{\sum_{\ell=1}^{J} \lambda_{\ell,n} \left[ \sum_{m=1}^{M} \hat{\varphi}_{\ell,n}^{(\alpha)}(Y_{m,n}) + (\alpha - 1)\kappa_n \right]^{\eta_n}} $$

$$ \Sigma_{j,n+1} = (1 - \gamma_n)\Sigma_{j,n} + \gamma_n \hat{\Sigma}_{j,n} + \gamma_n(1 - \gamma_n)(\hat{m}_{j,n} - m_{j,n})(\hat{m}_{j,n} - m_{j,n})^T. $$

**Remark 2** Approximated update rules for GMMs can be obtained by sampling according to $q_n = N(\cdot; 0_d, I_d)$, where $0_d$ is the null vector of dimension $d$. These can be deduced by applying the reparameterisation trick: for all $m = 1...M$, $Y_{m,n}^{(j)} \sim k(\theta_{j,n}, \cdot)$ if $Y_{m,n}^{(j)} = m_{j,n} + \Sigma_{j,n}^{-1/2} \varepsilon_{m,n}$ with $\varepsilon_{m,n} \sim q_n$ and by observing that an unbiased estimator of (77) is $M^{-1} \sum_{m=1}^{M} (\mu_n k(Y_{m,n}^{(j)})/p(Y_{m,n}^{(j)}))^{\alpha-1} g(Y_{m,n}^{(j)})$. This choice of $q_n$ relies on the existence of a reparameterisation and while it incurs a gain in terms of sampling cost as we only need M samples, the evaluation cost remains $M \times J$.

We thus have access to tractable algorithms that iteratively update both the weights and components parameters of a GMM by optimising the $\alpha$-divergence between the mixture distribution and the targeted distribution. Our framework also allows the use of learning rates $\gamma_{j,n}$ that are dependent on $j$. This means that the learning rates can differ for each component parameters update and it paves the way for adaptive learning rates per components. This aspect is left for future work as we will focus on the case $\gamma_{j,n} := \gamma_n$ according to Algorithm 3 and move on to presenting our numerical experiments.

**Remark 3** Computing the estimator (78) can be expensive for two reasons: (i) $\mu_n k$ involves a sum over $J$ so that evaluating this function requires heavier computations as the number of components $J$ grows and (ii) bayesian applications with $p = p(\cdot, \mathcal{D})$ often consider large-scale data sets $\mathcal{D}$. To alleviate this computational burden, one can (i) approximate the summation
appearing in $\mu_n k$ using subsampling and (ii) use a mini-batch approach to approximate $p = p(\cdot, \mathcal{D})$, as detailed in (Li and Turner, 2016, Section 4.3).

# 8. Numerical Experiments

In our numerical experiments, we are interested in understanding the behaviour of Algorithm 3 in practice. We first investigate the challenging case of multimodal targets, as our framework is designed to handle multimodality thanks to the hyperparameter $\alpha \in [0, 1)$ and to the use of mixture models as a variational family.

## 8.1 Toy multimodal targets

Let $u_d$ be the $d$-dimensional vector whose coordinates are all equal to 1, $I_d$ be the identity matrix and $c$ be a positive constant (we set $c = 2$). We consider three multimodal examples:

(i) **Equally-weighted Gaussian Mixture Model.** The target $p$ is a mixture density of two equally-weighted $d$-dimensional Gaussian distributions multiplied by $c$ such that

$$p(y) = c \times [0.5\mathcal{N}(y; -2u_d, I_d) + 0.5\mathcal{N}(y; 2u_d, I_d)] .$$

(ii) **Imbalanced Gaussian Mixture Model.** The target $p$ is a mixture density of three $d$-dimensional Gaussian distributions with unequal weights and multiplied by $c$ such that

$$p(y) = c \times [0.35\mathcal{N}(y; -2u_d, I_d) + 0.25\mathcal{N}(y; 2u_d, I_d) + 0.4\mathcal{N}(y; u_d, I_d)]$$

(iii) **Equally-weighted Student's t Mixture Model.** The target $p$ is a mixture density of two equally-weighted $d$-dimensional Student's $t$ distributions with two degree of freedom (i.e. $a = 2$) multiplied by $c$ such that

$$p(y) = c \times [0.5 t(y; -2u_d, I_d, a) + 0.5 t(y; 2u_d, I_d, a)] .$$

### 8.1.1 COMPARING THE RGD AND THE MG APPROACHES WITH $\eta_n = 0$

Setting $\eta_n = 0$ at all times $n \ge 1$ keeps the mixture weights fixed in Algorithm 3. In that case, our RGD approach relates to Gradient Ascent steps on the VR Bound (Li and Turner, 2016) w.r.t. the means of our Gaussian components (see Section 5.1). We then want to compare the performances of the RGD approach (that can be derived from Li and Turner, 2016) to the novel MG approach we have introduced in our work.

**Implementation details.** The covariance matrices of the mixture components are fixed and equal to $\sigma^2 I_d$ with $\sigma^2 = 1$, $\alpha = 0.2$, $J \in \{10, 50\}$, $d = 16$, $M = 200$, the total number of iterations $N$ is equal to 100, $\Theta_1$ is generated by sampling from a centered normal distribution with covariance matrix $10I_d$, $\lambda_1 = [1/J, \dots, 1/J]$ and for all time $n = 1 \dots N$, $\kappa_n = 0$, $\eta_n = 0$. and $\gamma_n = \gamma$ with $\gamma \in \{0.1, 0.5, 1.\}$. The experiments are replicated independently 30 times and the convergence of the RGD and of the MG approaches is monitored for the three multimodal examples (i), (ii) and (iii) by computing a Monte Carlo estimator of the VR Bound (since we have sampled $M$ samples from $q_n$ at time $n$, these samples can readily be used to obtain an estimate the VR Bound with no additional computational cost).
<figure><img src="image_5.png" /><figcaption>Figure 1: Monte Carlo estimate of the VR Bound for the RGD and the MG approaches (fixed mixture weights) when considering each of the target distributions (i), (ii) and (iii).</figcaption></figure>

Our results are plotted on Figure 1, where RGD-IS-n(γ) and MG-IS-n(γ) denote the algorithms originating from the RGD and MG approaches (the IS-n and IS-unif samplers are equivalent here), and error bounds plots can be found in Figures 5 and 6 of Appendix D. As minimising Ψ_α is equivalent to maximising the VR Bound when α ∈ (0, 1), our theoretical results predict a systematic increase in the VR Bound for the algorithms considered. This is what we observe in Figure 1 and we also see that the choice of γ plays the role of a learning rate our algorithms: while increasing γ mostly improves the speed of convergence, it may deteriorate the performances if chosen too large, such as when J = 50 with the target (iii).

Furthermore, MG-IS-n(γ) leads to better performances in all but one case compared to RDG-IS-n(γ). This renders the MG-IS-n(γ) algorithm more suitable overall for capturing the multimodality of the various multimodal targets p we considered, a fact that is further
<table><caption>Table 2: LogMSE averaged for the RGD and the MG approaches (fixed mixture weights) when considering each of the target distributions (i), (ii) and (iii).</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">J = 10</th><th colspan="3">J = 50</th></tr><tr><th>γ = 0.1</th><th>γ = 0.5</th><th>γ = 1.0</th><th>γ = 0.1</th><th>γ = 0.5</th><th>γ = 1.0</th></tr></thead><tbody><tr><td rowspan="2">(i)</td><td>RGD-IS-n(γ)</td><td>-0.081</td><td>-0.076</td><td>-0.218</td><td>-1.640</td><td>-1.673</td><td>-1.560</td></tr><tr><td>MG-IS-n(γ)</td><td><b>-3.702</b></td><td><b>-1.875</b></td><td><b>-2.711</b></td><td><b>-2.760</b></td><td><b>-2.771</b></td><td><b>-2.788</b></td></tr><tr><td rowspan="2">(ii)</td><td>RGD-IS-n(γ)</td><td>-0.211</td><td>-0.072</td><td>-0.015</td><td>-1.401</td><td>-1.437</td><td>-1.515</td></tr><tr><td>MG-IS-n(γ)</td><td><b>-2.581</b></td><td><b>-2.101</b></td><td><b>-1.742</b></td><td><b>-2.611</b></td><td><b>-2.328</b></td><td><b>-1.933</b></td></tr><tr><td rowspan="2">(iii)</td><td>RGD-IS-n(γ)</td><td>-0.108</td><td>-0.008</td><td>-0.111</td><td>-1.652</td><td>-1.654</td><td><b>-1.634</b></td></tr><tr><td>MG-IS-n(γ)</td><td><b>-0.913</b></td><td><b>-1.489</b></td><td><b>-1.846</b></td><td><b>-2.036</b></td><td><b>-2.530</b></td><td>-0.717</td></tr></tbody></table>

supported in Table 2, in which we compare the log of the Mean-Squared Error (logMSE) returned by each algorithm. (The MSE is computed as the average of $||m_{approx} - m_{true}||^2$ over 30 independent runs of each algorithm, where $||.||$ stands for the Euclidean norm, $m_{true} = E_{\mathbb{P}[Y]}$ with $\mathbb{P}(A) = \int_A p(y)\nu(dy) / \int_Y p(y)\nu(dy)$ for all $A \in \mathcal{Y}$ and where $m_{approx}$ is the mean of the optimised variational distribution, that is in our setting $m_{approx} = \sum_{j=1}^J \lambda_{j,N} m_{j,N}$.)

Let us next pair up the means optimisation with mixture weights optimisation and investigate the impact of the sequence $(\eta_n)_{n \ge 1}$ on our algorithms.

### 8.1.2 COMPARING THE RGD AND THE MG APPROACHES WITH $\eta_n > 0$.

So far, we have demonstrated the viability of our novel MG approach compared to the more usual Rényi's $\alpha$-divergence-based RGD approach. Yet, we have not taken advantage of the fact that we can perform mixture weights optimisation on top of means optimisation, which is another novelty of our framework compared to traditional Variational Inference methods.

*Implementation details.* We work with the same implementation details as those described in Section 8.1.1, except for the fact that we will now let $\eta_n = \eta$ at time $n$ and we will vary the value of $\eta$. In addition to the RGD-IS-n($\gamma$) and the MG-IS-n($\gamma$) algorithms, we also include the RGD-IS-unif($\gamma$) and the MG-IS-unif($\gamma$) algorithms in our results (that is, the algorithms resulting from pairing up the RGD and MG approaches with the IS-unif sampler).

The plots for the averaged Monte Carlo estimate of the VR Bound with $\eta = 0.1$ and $\gamma = 0.5$ are then available in Figure 2 (and similar plots for $\eta \in \{0.05, 0.5\}$ are provided in Figures 7 and 8 of Appendix D respectively). According to Figure 2, MG-IS-unif(0.5) outperforms most of the time the other methods (this trend for MG-IS-unif($\gamma$) is further confirmed by looking at the logMSE in Table 3 of Appendix D, which includes the cases $\gamma \in \{0.1, 1.\}$ for completeness).

We also plot in Figure 3 the final mixture weights $\lambda_N$ obtained after one typical run of MG-IS-unif(0.5) for $\eta \in \{0.05, 0.1, 0.5\}$ and $J \in \{10, 50\}$ (and additional LogMSE results can be found in Table 4 of Appendix D). This is done in order to delve into how $\eta$ impacts the variational density returned at the end of Algorithm 3 (MG-IS-unif(0.5) was chosen among all four options since it enjoys good empirical performances according to Figure 2).
<figure><img src="image_6.png" /><figcaption>Figure 2: Monte Carlo estimate of the VR Bound for the RGD and the MG approaches ($\eta = 0.1$) when considering each of the target distributions (i), (ii) and (iii).</figcaption></figure>

A key insight from Figure 3 is that optimising the mixture weights permits us to select the components appearing in our mixture model according to their overall contribution in constructing a good approximation of the targeted density.

On the one hand, we enable more flexibility in our variational approximation as we optimise over a set of component parameters $\Theta$ with $J > 1$ instead of just one component parameter $\theta$ ($J = 1$). On the other hand, we avoid complexifying unnecessarily the variational distribution returned at the end of the optimisation procedure since we perform mixture weights optimisation alongside components parameters optimisation. That way, we can bypass the limitation of the case $\eta = 0$, which may use more components than needed when approximating the targeted density (this is the case here since the targeted densities considered have three modes at best while $J = \{10, 50\}$).

However, there is a tradeoff to find between the simplicity of the variational density that is returned and its accuracy at describing the targeted density, which is expressed via
<figure><img src="image_7.png" /><figcaption>Figure 3: Final mixture weights $\lambda_N$ for one run of the MG-IS-unif(0.5) algorithm (varying mixture weights) when considering each of the target distributions (i), (ii) and (iii).</figcaption></figure>

the choice of the learning rate $\eta$. Indeed, choosing $\eta$ too large might lead to missing some of the modes while having $\eta$ too small may not discriminate quickly enough between the components. Observe in particular that the number of components $J$ too plays a role in how fast the components selection process occurs, as for the same value of $\eta$, this selection process is slower as $J$ increases.

## 8.2 Bayesian Logistic Regression

We consider the Bayesian Logistic Regression setting from Gershman et al. (2012). Namely, the data $\mathcal{D} = \{c, x\}$ is made of $I$ binary class labels, $c_i \in \{-1, 1\}$, and of $L$ covariates for each datapoint, $x_i \in \mathbb{R}^L$. The latent variable $y = \{w, \beta\}$ consists of $L$ regression coefficients
<figure><img src="image_8.png" /><figcaption>Figure 4: Monte Carlo estimate of the VR Bound for the RGD and the MG approaches when considering the Bayesian Logistic Regression on the Covertype data set.</figcaption></figure>

$w_l \in \mathbb{R}$, and a precision parameter $\beta \in \mathbb{R}^+$. Furthermore, the following model is assumed

$$p_0(\beta) = \text{Gamma}(\beta; a, b), \quad p_0(\mathbf{w}|\beta) = \mathcal{N}(\mathbf{w}; 0, \beta^{-1}\mathbf{I}_L), \quad p(c_i = 1|\mathbf{x}_i, \mathbf{w}) = \frac{1}{1 + e^{-\mathbf{w}^T\mathbf{x}_i}},$$

where $a = 1$ and $b = 0.01$, so that $p(y, \mathcal{D}) = p_0(\beta)p_0(\mathbf{w}|\beta) \prod_{i=1}^{I} p(c_i|\mathbf{x}_i, y)$. We now want to demonstrate the practicability of our framework in a real data scenario. To this end, we select the Covertype data set (581,012 data points and 54 features, available at https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html) and we compare the RGD and the MG approaches for this choice of target density.

*Implementation details.* The covariance matrices of the mixture components are fixed and equal to $\sigma^2 \mathbf{I}_d$ with $\sigma^2 = 1$, $\alpha = 0.2$, $J \in \{50, 100\}$, $d = 56$, $M = 200$, the total number of iterations $N$ is equal to 200, $\Theta_1$ is generated by sampling from a centered normal distribution with covariance matrix $5\mathbf{I}_d$, $\lambda_1 = [1/J, \dots, 1/J]$ and for all time $n = 1 \dots N$, $\kappa_n = 0$, $\eta_n = \eta$, $\gamma_n = \gamma$ where $\eta \in \{0.1, 0.5\}$ and $\gamma \in \{0.05, 0.1, 0.2, 0.3\}$. Computing $p(y, \mathcal{D})$ constitutes the major computation bottleneck here since $I = 581,012$. To address this problem, $p(y, \mathcal{D})$ is approximated according to Remark 3 with subsampled mini-batches of size 100. The experiments are replicated independently 30 times and the convergence of the RGD and of the MG approaches is monitored by computing a Monte Carlo estimator of the VR Bound.

Our results are plotted in Figure 4, in which we focus on RGD-IS-n($\gamma$) and MG-IS-unif($\gamma$) (those two versions of Algorithm 3 are the most interesting to compare in this particular setting since MG-IS-n($\gamma$) enjoys similar performances to RGD-IS-n($\gamma$) and RGD-IS-unif($\gamma$)
underperforms, see Figure 9 of Appendix D for details). We observe that both algorithms are able to learn in a real data scenario for a proper tuning of $(\eta, \gamma)$ and that selecting either $\eta$ or $\gamma$ too small/large can deteriorate the performance (as per the learning rate behaviour associated to those parameters). Furthermore, MG-IS-unif($\gamma$) can be tuned to outperform RGD-IS-n($\gamma$), which confirms our previous empirical findings regarding the relevance of the novel MG approach compared to the more traditional RGD one.

## 9. Conclusion

We introduced a novel methodology to build algorithms ensuring a monotonic decrease in the $\alpha$-divergence at each step. Our methodology enabled simultaneous updates for both the weights and components parameters of a given mixture model, making it suitable for capturing complex multimodal target densities. Our work also connected and improved on different approaches: Gradient Descent for $\alpha$-divergence minimisation, Power Descent and an Integrated EM algorithm. By investigating variational families based on the exponential family, we applied our framework to important classes of models such as Gaussian Mixture Models and Student's $t$ Mixture Models. Finally, we provided empirical evidence that our methodology can be used to enhance the aforementioned existing algorithms.

To conclude, we state several directions to extend our work. Now that we have established a systematic decrease for our iterative schemes, the next step could be to derive convergence results and to compare them with those obtained using typical Gradient Descent schemes. Based on Proposition 2, another direction is to generalise the monotonicity property from Theorem 3 beyond the case $\alpha \in [0, 1)$. Lastly, much like it is already the case in traditional Black-Box Variational Inference, we expect that fine-tuning our hyperparameters and resorting to more advanced Monte Carlo methods in the estimation of the intractable integrals appearing in our ideal algorithms will lead to further improved numerical results.

## Acknowledgments and Disclosure of Funding

We would like to thank the action editor and the reviewers for helpful comments and suggestions on the paper. Kamélia Daudel acknowledges support of the UK Defence Science and Technology Laboratory (Dstl) and and Engineering and Physical Research Council (EPSRC) under grant EP/R013616/1. This is part of the collaboration between US DOD, UK MOD and UK EPSRC under the Multidisciplinary University Research Initiative.

## A. Deferred Proofs and Results for Section 3

### A.1 Quantifying the Improvement in one Step of Gradient Descent

Let $T \subseteq \mathbb{R}^d$ be a convex set. Here, $\langle \cdot, \cdot \rangle$ is the standard inner product on $\mathbb{R}^d$ and $\|\cdot\|$ is the Euclidean norm.

**Definition 3** A continuously differentiable function $g$ defined on $T$ is said to be $\beta$-smooth if for all $\theta, \theta' \in T$,

$$ \|\nabla g(\theta) - \nabla g(\theta')\| \le \beta \|\theta - \theta'\|. $$
**Lemma 2** Let $\gamma \in (0, 1]$, let $g$ be a $\beta$-smooth function defined on $\mathcal{T}$. Then, for all $\theta \in \mathcal{T}$ it
holds that

$$
g(\theta) - g\left(\theta - \frac{\gamma}{\beta} \nabla g(\theta)\right) \geq \frac{\gamma}{2\beta} \|\nabla g(\theta)\|^2.
$$

**Proof** By assumption on $g$, we have that for all $\theta, \theta' \in \mathcal{T}$

$$
g(\theta') - g(\theta) - \langle \nabla g(\theta), \theta' - \theta \rangle \le \frac{\beta}{2} \|\theta - \theta'\|^2.
$$

In particular, setting $\theta' = \theta - \frac{\gamma}{\beta} \nabla g(\theta)$ yields

$$
\begin{align*}
g(\theta) - g\left(\theta - \frac{\gamma}{\beta} \nabla g(\theta)\right) &\ge \frac{\gamma}{\beta} \|\nabla g(\theta)\|^2 - \frac{\gamma^2}{2\beta} \|\nabla g(\theta)\|^2 \\
&\ge \frac{\gamma}{\beta} \left(1 - \frac{\gamma}{2}\right) \|\nabla g(\theta)\|^2.
\end{align*}
$$

Since $\gamma \in (0, 1]$, we can deduce the desired result.

<figure><img src="image_9.png" /></figure>

**A.2 Gradient Descent for α- / Rényi's α-divergence Minimisation in Section 3**

Let us first write the definition of the α-divergence (resp. of Rényi's α-divergence) between
the two absolutely continuous probability measures $K(\theta, \cdot)$ and $\mathbb{P}$

$$
\begin{align*}
D_{\alpha}(K(\theta, \cdot) || \mathbb{P}) &= \int_{\mathcal{Y}} \frac{1}{\alpha(\alpha - 1)} \left[ \left( \frac{k(\theta, y)}{p(y|\mathcal{D})} \right)^{\alpha} - 1 \right] p(y|\mathcal{D}) \nu(dy) \\
D_{\alpha}^{(\text{AR})}(K(\theta, \cdot) || \mathbb{P}) &= \frac{1}{\alpha(\alpha - 1)} \log \left( \int_{\mathcal{Y}} k(\theta, y)^{\alpha} p(y|\mathcal{D})^{1-\alpha} \nu(dy) \right).
\end{align*}
$$

Here, the Rényi divergence is defined following the convention from Cichocki and Amari
(2010), alternative definitions may use a different scaling factor.

A.2.1 GRADIENT DESCENT FOR $\alpha$-DIVERGENCE MINIMISATION

As underlined in the introduction, minimising the α-divergence $D_\alpha(\mathbb{Q}||\mathbb{P})$ w.r.t $q$ is equivalent
to minimising $\Psi_\alpha(q;p)$ with $p = p(\cdot, \mathcal{D})$ w.r.t $q$, where we have gotten rid of $p(\mathcal{D})$ in the opti-
misation problem as written in (2). Letting $\mathcal{Q}$ be of the form $\mathcal{Q} = \{q : y \mapsto k(\theta, y) : \theta \in \mathcal{T}\}$,
the traditional Variational Inference way to optimise $\Psi_\alpha(k(\theta, \cdot); p)$ with $p = p(\cdot, \mathcal{D})$ w.r.t
$\theta$ corresponds to performing Gradient Descent steps on $\theta$ to construct a sequence $(\theta_n)_{n \ge 1}$
that converges towards a local optimum of the function $\theta \mapsto \Psi_\alpha(k(\theta, \cdot); p)$. This procedure
involves a well-chosen learning rate policy $(r_n)_{n \ge 1}$ and sets:

$$
\theta_{n+1} = \theta_n - r_n \nabla \Psi_{\alpha}(k(\theta_n, \cdot); p), \quad n \ge 1. \tag{79}
$$
Under common differentiability assumptions,

$$
\begin{align*}
\nabla \Psi_{\alpha}(k(\theta_n, \cdot); p) &= \nabla \left( \int_{\mathcal{Y}} f_{\alpha} \left( \frac{k(\theta_n, y)}{p(y)} \right) p(y) \nu(dy) \right) \\
&= \int_{\mathcal{Y}} \left. \frac{\partial}{\partial \theta} \left( f_{\alpha} \left( \frac{k(\theta, y)}{p(y)} \right) \right) \right|_{(\theta, y)=(\theta_n, y)} p(y) \nu(dy) \\
&= \int_{\mathcal{Y}} f'_{\alpha} \left( \frac{k(\theta_n, y)}{p(y)} \right) \left. \frac{\partial k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy) \\
&= \int_{\mathcal{Y}} \frac{k(\theta_n, y)^{\alpha} p(y)^{1-\alpha}}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy)
\end{align*}
$$

so that (79) becomes

$$
\theta_{n+1} = \theta_n - r_n \int_{\mathcal{Y}} \frac{\varphi_n^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy), \quad n \ge 1.
$$

A.2.2 GRADIENT DESCENT FOR RÉNYI’S α-DIVERGENCE MINIMISATION

Considering yet again $\mathcal{Q} = \{q : y \mapsto k(\theta, y) : \theta \in \mathcal{T}\}$, minimising Rényi's $\alpha$-divergence

$$
D_{\alpha}^{(\mathrm{AR})}(K(\theta, \cdot)\|\mathbb{P}) = \frac{1}{\alpha(\alpha - 1)} \log \left( \int_{\mathcal{Y}} k(\theta, y)^{\alpha} p(y|\mathcal{D})^{1-\alpha} \nu(dy) \right)
$$

w.r.t $\theta$ can be done by performing Gradient Descent steps on $\theta$ to construct a sequence $(\theta_n)_{n \ge 1}$ that converges towards a local optimum of the function $\theta \mapsto D_{\alpha}^{(\text{AR})}(K(\theta, \cdot)\|\mathbb{P})$. This procedure involves a well-chosen learning rate policy $(r_n)_{n \ge 1}$ and sets:

$$
\theta_{n+1} = \theta_n - r_n \nabla D_{\alpha}^{(\text{AR})}(K(\theta_n, \cdot) || \mathbb{P}), \quad n \ge 1. \tag{80}
$$

Under common differentiability assumptions and setting $p = p(\cdot, \mathcal{D})$,

$$
\begin{align*}
\nabla D_{\alpha}^{(\text{AR})}(K(\theta_n, \cdot) || \mathbb{P}) &= \nabla \left( \frac{1}{\alpha(\alpha-1)} \log \left( \int_{\mathcal{Y}} k(\theta_n, y)^{\alpha} p(y|\mathcal{D})^{1-\alpha} \nu(dy) \right) \right) \\
&= \frac{1}{\alpha(\alpha-1)} \frac{\nabla \left( \int_{\mathcal{Y}} k(\theta_n, y)^{\alpha} p(y, \mathcal{D})^{1-\alpha} \nu(dy) \right)}{\int_{\mathcal{Y}} k(\theta_n, y)^{\alpha} p(y, \mathcal{D})^{1-\alpha} \nu(dy)} \\
&= \int_{\mathcal{Y}} \frac{\tilde{\varphi}_n^{(\alpha)}(y)}{\alpha-1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy)
\end{align*}
$$

so that (80) becomes

$$
\theta_{n+1} = \theta_n - r_n \int_{\mathcal{Y}} \frac{\tilde{\varphi}_n^{(\alpha)}(y)}{\alpha - 1} \left. \frac{\partial \log k(\theta, y)}{\partial \theta} \right|_{(\theta, y)=(\theta_n, y)} \nu(dy), \quad n \ge 1.
$$

B. Deferred Proofs and Results for Section 4

B.1 Proof of Corollary 4

**Proof of Corollary 4** Following the proof of Corollary 1, we obtain that (19) holds by using the definition of $\Theta_{n+1}$ combined with the fact that $\alpha \in [0, 1)$, $D_1(K(\theta_{j,n}, \cdot) || K(\theta_{j,n+1}, \cdot)) \ge 0$
and $\lambda_{j,n} > 0$ for all $j = 1...J$. (18) holds by Theorem 3 and we can thus apply Theorem 2. ■

## B.2 Extension of Lemma 1 to Mixture Models

**Lemma 3 (Generalised maximisation approach for the mean-field family)** Let each component of the mixture be a member of the same mean-field variational family so that $k(\theta_j, y) = \prod_{\ell=1}^{L} k^{(\ell)}(\theta_j^{(\ell)}, y^{(\ell)})$ with $\theta_j = (\theta_j^{(1)},...,\theta_j^{(L)}) \in \mathcal{T}$. Then, starting from $\Theta_1 \in \mathcal{T}$ and denoting $\Theta_n = (\theta_{1,n},...,\theta_{j,n})$ with $\theta_{j,n} = (\theta_{j,n}^{(1)},...,\theta_{j,n}^{(L)})$ for all $n \ge 1$ and all $j = 1...J$, solving (25) yields the following update formulas: for all $n \ge 1$ and all $j = 1...J$,

$$ \theta_{j,n+1}^{(\ell)} = \underset{\theta^{(\ell)}}{\operatorname{argmax}} \int_{\mathcal{Y}} [\varphi_{j,n}^{(\alpha)}(y) + b_{j,n}k(\theta_n, y)] \log \left( \frac{k^{(\ell)}(\theta^{(\ell)}, y^{(\ell)})}{k^{(\ell)}(\theta_{j,n}^{(\ell)}, y^{(\ell)})} \right) \nu(dy), \quad \ell = 1...L. $$

## B.3 Proof of Corollary 5

**Proof of Corollary 5** Following the proof of Corollary 2, we will use that $\gamma_{j,n} \in (0,1]$ and that $g_{j,n}$ is a $\beta_{j,n}$-smooth function. Indeed, we can thus apply Lemma 2 and we obtain by definition of $\theta_{j,n+1}$ in (26) that $0 = g_{j,n}(\theta_{j,n}) \ge g_{j,n}(\theta_{j,n+1})$ for all $n \ge 1$ and all $j = 1...J$, which in turn implies (19). In addition, (18) holds by Theorem 3, hence we can apply Theorem 2. ■

## B.4 Gradient Descent for $\alpha$- / Rényi's $\alpha$-divergence Minimisation in Section 4

We obtain the desired results by adapting the reasoning from Appendix A.2 to the case

$$ Q = \left\{ q : y \mapsto \mu_{\lambda, \Theta} k(y) = \sum_{j=1}^{J} \lambda_j k(\theta_j, y) : \Theta \in \mathcal{T}^{J} \right\}, $$

with $\lambda \in S_J$, meaning that we consider the $\alpha$-divergence (resp. Rényi's $\alpha$-divergence) between the two absolutely continuous probability measures $\mu_{\lambda, \Theta} K$ and $\mathbb{P}$:

$$ D_{\alpha}(\mu_{\lambda, \Theta} K || \mathbb{P}) = \int_{\mathcal{Y}} \frac{1}{\alpha(\alpha - 1)} \left[ \left( \frac{\mu_{\lambda, \Theta} k(y)}{p(y|\mathcal{D})} \right)^{\alpha} - 1 \right] p(y|\mathcal{D}) \nu(dy) $$

$$ D_{\alpha}^{(\text{AR})}(\mu_{\lambda, \Theta} K || \mathbb{P}) = \frac{1}{\alpha(\alpha - 1)} \log \left( \int_{\mathcal{Y}} (\mu_{\lambda, \Theta} k(y))^{\alpha} p(y|\mathcal{D})^{1-\alpha} \nu(dy) \right), $$

(we use the convention from Cichocki and Amari (2010) for Rényi's $\alpha$-divergence).

## B.5 Monotonicity Property for the Power Descent

### B.5.1 PRELIMINARY REMARKS

For convenience, we redefine in this section the function $b_{\mu,\alpha}$ for all $\mu \in M_1(\mathcal{T})$ by

$$ b_{\mu, \alpha}(\theta) = \int_{\mathcal{Y}} k(\theta, y) \frac{1}{\alpha - 1} \left( \frac{\mu k(y)}{p(y)} \right)^{\alpha - 1} \nu(dy), \quad \theta \in \mathcal{T}. $$
Then, for all $\eta > 0$, the iteration $\mu \mapsto \mathcal{I}_\alpha(\mu)$ is well-defined if we have

$$0 < \mu(|b_{\mu,\alpha} + \kappa|^{\frac{\eta}{1-\alpha}}) < \infty. \qquad (81)$$

Furthermore, Daudel et al. (2021) already established that one transition of the Power Descent algorithm ensures a monotonic decrease in the $\alpha$-divergence at each step for all $\eta \in (0, 1]$ and all $\kappa$ such that $(\alpha - 1)\kappa \ge 0$ under the assumption of Proposition 2, which settles the case (iii). Finally, while we establish our results for (i) and (ii) in the general case where $\mu \in M_1(\mathcal{T})$, the particular case of mixture models follows immediately by choosing $\mu$ as a weighted sum of dirac measures.

### B.5.2 EXTENDING THE MONOTONICITY

Let $(\zeta, \mu)$ be a couple of probability measures where $\zeta$ is dominated by $\mu$, which we denote by $\zeta \preceq \mu$. A first lower-bound for the difference $\Psi_\alpha(\mu k) - \Psi_\alpha(\zeta k)$ was derived in Daudel et al. (2021) and was used to establish that the Power Descent algorithm diminishes $\Psi_\alpha$ for all $\eta \in (0, 1]$. We now prove a novel lower-bound for the difference $\Psi_\alpha(\mu k) - \Psi_\alpha(\zeta k)$ which will allow us to extend the monotonicity results from Daudel et al. (2021) beyond the case $\eta \in (0, 1]$ when $\alpha < 0$. This result relies on the existence of an exponent $\varrho$ satisfying condition (A2) below, which will later on be used to specify a range of values for $\eta$ ensuring that $\Psi_\alpha$ is decreasing after having applied one transition $\mu \mapsto \mathcal{I}_\alpha(\mu)$

(A2) We have $\varrho \in \mathbb{R} \setminus [0, 1]$ and the function $f_{\alpha, \varrho} : u \mapsto f_\alpha(u^{1/\varrho})$ is non-decreasing and concave on $\mathbb{R}_{>0}$.

**Proposition 5** Assume (A1). Let $\alpha \in \mathbb{R} \setminus \{1\}$, assume that $\varrho$ satisfies (A2) and let $\kappa$ be such that $(\alpha - 1)\kappa \ge 0$. Then, for all $\mu, \zeta \in M_1(\mathcal{T})$ such that $\mu(|b_{\mu,\alpha}|) < \infty$ and $\zeta \preceq \mu$,

$$|\varrho|^{-1} \{\mu(|b_{\mu,\alpha} + \kappa|) - \mu(|b_{\mu,\alpha} + \kappa|g^{\varrho})\} \le \Psi_{\alpha}(\mu k) - \Psi_{\alpha}(\zeta k), \qquad (82)$$

where $g$ is the density of $\zeta$ wrt $\mu$, i.e. $\zeta(d\theta) = \mu(d\theta)g(\theta)$. Moreover, equality holds in (82) if and only if $\zeta = \mu$.

**Proof** First note that for all $\alpha \in \mathbb{R} \setminus \{1\}$, we have by (A2) that $f'_{\alpha, \varrho}(u) \ge 0$ for all $u > 0$, and thus that $sg(\varrho) = sg(\alpha - 1)$ where $sg(v) = 1$ if $v \ge 0$ and $-1$ otherwise. Since $sg(f'_\alpha(u)) = sg(\alpha - 1) = sg(\kappa)$ for all $u > 0$, this implies that $\varrho^{-1}f'_\alpha(u) = |\varrho|^{-1}|f'_\alpha(u)|$, $\varrho^{-1}\kappa = |\varrho^{-1}\kappa|$ and finally that $\varrho^{-1}(b_{\mu,\alpha}(\theta) + \kappa) = |\varrho^{-1}||b_{\mu,\alpha}(\theta) + \kappa|$ for all $\theta \in \mathcal{T}$, which will be used later in the proof.

Write by definition of $f_{\alpha, \varrho}$ in (A2) and $\zeta$,

$$ \begin{aligned} \Psi_{\alpha}(\zeta k) &= \int_{\mathcal{Y}} f_{\alpha} \left( \frac{\zeta k(y)}{p(y)} \right) p(y) \nu(dy) \\ &= \int_{\mathcal{Y}} f_{\alpha, \varrho} \left( \left[ \frac{\zeta k(y)}{p(y)} \right]^{\varrho} \right) p(y) \nu(dy) \\ &= \int_{\mathcal{Y}} f_{\alpha, \varrho} \left( \left[ \int_{\mathcal{T}} \mu(d\theta) \frac{k(\theta, y)}{\mu k(y)} \left( \frac{g(\theta) \mu k(y)}{p(y)} \right) \right]^{\varrho} \right) p(y) \nu(dy) \\ &\le \int_{\mathcal{Y}} f_{\alpha, \varrho} \left( \int_{\mathcal{T}} \mu(d\theta) \frac{k(\theta, y)}{\mu k(y)} \left( \frac{g(\theta) \mu k(y)}{p(y)} \right)^{\varrho} \right) p(y) \nu(dy) \end{aligned} \qquad (83) $$
where the last inequality follows from Jensen's inequality applied to the convex function $u \mapsto u^\varrho$ (since $\varrho \in \mathbb{R} \setminus [0, 1]$) and the fact that $f_{\alpha, \varrho}$ is non-decreasing. Now set

$$ 
\begin{aligned}
u_y &= \int_{\mathbb{T}} \mu(d\theta) \frac{k(\theta, y)}{\mu k(y)} \left( \frac{g(\theta) \mu k(y)}{p(y)} \right)^{\varrho} \\
v_y &= \left( \frac{\mu k(y)}{p(y)} \right)^{\varrho}
\end{aligned}
$$

and note that

$$ u_y - v_y = \left( \frac{\mu k(y)}{p(y)} \right)^{\varrho} \left( \int_{\mathbb{T}} \mu(d\theta) \frac{k(\theta, y)}{\mu k(y)} g^{\varrho}(\theta) - 1 \right). \quad (84) $$

Since $f_{\alpha, \varrho}$ is concave, $f_{\alpha, \varrho}(u_y) \le f_{\alpha, \varrho}(v_y) + f'_{\alpha, \varrho}(v_y)(u_y - v_y)$. Combining with (83), we get

$$ 
\begin{aligned}
\Psi_{\alpha}(\zeta k) &\le \int_{\mathbb{Y}} f_{\alpha, \varrho}(u_y)p(y)\nu(dy) \\
&\le \int_{\mathbb{Y}} f_{\alpha, \varrho}(v_y)p(y)\nu(dy) + \int_{\mathbb{Y}} f'_{\alpha, \varrho}(v_y)(u_y - v_y)p(y)\nu(dy)
\end{aligned}
\quad (85) $$

Note that the first term of the rhs can be written as

$$ \int_{\mathbb{Y}} f_{\alpha, \varrho}(v_y)p(y)\nu(dy) = \int_{\mathbb{Y}} f_{\alpha}\left(\frac{\mu k(y)}{p(y)}\right)p(y)\nu(dy) = \Psi_{\alpha}(\mu k). \quad (86) $$

Using now $f'_{\alpha, \varrho}(v_y) = \varrho^{-1}v_y^{1/\varrho-1}f'_{\alpha}(v_y^{1/\varrho})$ and (84), the second term of the rhs of (85) may be expressed as

$$ 
\begin{aligned}
&\int_{\mathbb{Y}} f'_{\alpha, \varrho}(v_y)(u_y - v_y)p(y)\nu(dy) \\
&= \varrho^{-1} \int_{\mathbb{Y}} \left(\frac{\mu k(y)}{p(y)}\right)^{1-\varrho} f'_{\alpha}\left(\frac{\mu k(y)}{p(y)}\right) \\
&\qquad \left(\frac{\mu k(y)}{p(y)}\right)^{\varrho} \left(\int_{\mathbb{T}} \mu(d\theta) \frac{k(\theta, y)}{\mu k(y)} g^{\varrho}(\theta) - 1\right) p(y)\nu(dy) \\
&= \varrho^{-1} \int_{\mathbb{T}} \mu(d\theta) \left(\int_{\mathbb{Y}} k(\theta, y) f'_{\alpha}\left(\frac{\mu k(y)}{p(y)}\right) \nu(dy)\right) g^{\varrho}(\theta) \\
&\qquad - \varrho^{-1} \int_{\mathbb{Y}} \mu k(y) f'_{\alpha}\left(\frac{\mu k(y)}{p(y)}\right) \nu(dy) \\
&= \varrho^{-1} \left\{\mu(b_{\mu,\alpha} \cdot g^{\varrho}) - \mu(b_{\mu,\alpha})\right\} \\
&= |\varrho|^{-1} \left\{\mu(|b_{\mu,\alpha} + \kappa|g^{\varrho}) - \mu(|b_{\mu,\alpha} + \kappa|)\right\} + |\varrho|^{-1}\kappa(1 - \mu(g^{\varrho})),
\end{aligned}
$$

where we have used that $\varrho^{-1}(b_{\mu,\alpha}(\theta) + \kappa) = |\varrho^{-1}||b_{\mu,\alpha}(\theta) + \kappa|$ for all $\theta \in \mathbb{T}$ and that $\varrho^{-1}\kappa = |\varrho^{-1}\kappa|$. In addition, Jensen's inequality applied to the convex function $u \mapsto u^\varrho$ implies that $\mu(g^{\varrho}) \ge 1$ and thus

$$ \int_{\mathbb{Y}} f'_{\alpha, \varrho}(v_y)(u_y - v_y)p(y)\nu(dy) \le |\varrho|^{-1} \{\mu(|b_{\mu,\alpha} + \kappa|g^{\varrho}) - \mu(|b_{\mu,\alpha} + \kappa|)\}. \quad (87) $$
Combining this inequality with (85) and (86) finishes the proof of the inequality. Furthermore, if the equality holds in (82), then the equality in Jensen's inequality (87) shows that $g$ is constant $\mu$-a.e. so that $\zeta = \mu$, and the proof is completed. ■

**Remark 4** The proof of Proposition 5 relies on $f'_\alpha$ being of constant sign. Notice however that the definition of the $\alpha$-divergence in (1) is invariant with respect to the transformation $\tilde{f}_\alpha(u) = f_\alpha(u) + \kappa(u-1)$ for any arbitrary constant $\kappa$, that is $f_\alpha$ can be equivalently replaced by $\tilde{f}_\alpha$ in (1). This aspect is in fact expressed through the constant $\kappa$ appearing in the update formula, that we however need to assume to satisfy $(\alpha - 1)\kappa \ge 0$ in our proofs.

We now plan on setting $\zeta = \mathcal{I}_\alpha(\mu)$ in Proposition 5 and obtain that one iteration of the Power Descent yields $\Psi_\alpha(\mathcal{I}_\alpha(\mu)k) \le \Psi_\alpha(\mu k)$. For this purpose and based on the upper bound obtained in Proposition 5, we strengthen the condition (81) as follows to take into account the exponent $\rho$

$$0 < \mu(|b_{\mu,\alpha} + \kappa|^{\frac{\eta}{1-\alpha}}) < \infty \quad \text{and} \quad \mu(|b_{\mu,\alpha} + \kappa|^{\rho}) \le \mu(|b_{\mu,\alpha} + \kappa|) \\ \text{with } g = \frac{|b_{\mu,\alpha} + \kappa|^{\frac{\eta}{1-\alpha}}}{\mu(|b_{\mu,\alpha} + \kappa|^{\frac{\eta}{1-\alpha}})} \quad (88)$$

This leads to the following result.

**Proposition 6** Assume (A1). Let $\alpha \in \mathbb{R} \setminus \{1\}$, assume that $\rho$ satisfies (A2) and let $\kappa$ be such that $(\alpha - 1)\kappa \ge 0$. Let $\mu \in M_1(\mathbb{T})$ be such that $\mu(|b_{\mu,\alpha}|) < \infty$ and assume that $\eta$ satisfies (88). Then, the two following assertions hold.

(i) We have $\Psi_\alpha(\mathcal{I}_\alpha(\mu)k) \le \Psi_\alpha(\mu k)$.

(ii) We have $\Psi_\alpha(\mathcal{I}_\alpha(\mu)k) = \Psi_\alpha(\mu k)$ if and only if $\mu = \mathcal{I}_\alpha(\mu)$.

**Proof** We treat the case $\kappa = 0$ in the proof below (the case $\kappa \ne 0$ unfolds similarly). We apply Proposition 5 with $\zeta = \mathcal{I}_\alpha(\mu)$ so that $\zeta(d\theta) = \mu(d\theta)g(\theta)$ with $g = |b_{\mu,\alpha}|^{\eta/(1-\alpha)}/\mu(|b_{\mu,\alpha}|^{\eta/(1-\alpha)})$. Then,

$$\Psi_{\alpha}(\mathcal{I}_{\alpha}(\mu)k) \le \Psi_{\alpha}(\mu k) + |\rho|^{-1} \{\mu(|b_{\mu,\alpha}|g^{\rho}) - \mu(|b_{\mu,\alpha}|)\} \le \Psi_{\alpha}(\mu k), \quad (89)$$

where the last inequality follows from condition (88). Let us now show (ii). The *if* part is obvious. As for the *only if* part, $\Psi_\alpha(\mathcal{I}_\alpha(\mu)k) = \Psi_\alpha(\mu k)$ combined with (89) yields

$$\Psi_{\alpha}(\mathcal{I}_{\alpha}(\mu)k) = \Psi_{\alpha}(\mu k) + |\rho|^{-1} \{\mu(|b_{\mu,\alpha}|g^{\rho}) - \mu(|b_{\mu,\alpha}|)\},$$

which is the case of equality in Proposition 5. Therefore, $\mathcal{I}_\alpha(\mu) = \mu$. ■

While Proposition 6 resembles (Daudel et al., 2021, Theorem 1) in its formulation and in the properties on the iteration $\mu \mapsto \mathcal{I}_\alpha(\mu)$ it establishes, it is important to note that the proof techniques used, and thus the conditions on $\eta$ obtained, are different. This brings us to the proof of Proposition 2. The proof of this theorem requires intermediate results, which are derived in Appendix B.5.3 alongside the proof of Proposition 2.
### B.5.3 PROOF OF PROPOSITION 2

For the sake of readability, we only treat the case $κ = 0$ in the proofs below (the case $κ \neq 0$ unfolds similarly). In Proposition 5, the difference $\Psi_{\alpha}(\zeta k) - \Psi_{\alpha}(\mu k)$ is split into two terms

$$ \Psi_{\alpha}(\zeta k) - \Psi_{\alpha}(\mu k) = A(\mu, \zeta) + |\varrho|^{-1} \{\mu(|b_{\mu,\alpha}|g^{\varrho}) - \mu(|b_{\mu,\alpha}|)\}, $$

where $g = d\zeta/d\mu$. Moreover, Proposition 5 states that $A(\mu, \zeta)$ 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_{\mu,\alpha}|^{1/(1-\varrho)}$, as we show in Lemma 4 below.

**Lemma 4** Let $\varrho \in \mathbb{R} \setminus [0, 1]$. Then, for any positive probability density $g$ w.r.t $\mu$, we have

$$ \mu(|b_{\mu,\alpha}|g^{\varrho}) \geq [\mu(|b_{\mu,\alpha}|^{1/(1-\varrho)})]^{1-\varrho}, $$

with equality if and only if $g \propto |b_{\mu,\alpha}|^{1/(1-\varrho)}$.

**Proof** The function $x \mapsto x^{1-\varrho}$ is strictly convex for $\varrho \in \mathbb{R} \setminus [0, 1]$. Thus Jensen's inequality yields, for any positive probability density $g$ w.r.t. $\mu$,

$$ \mu(|b_{\mu,\alpha}|g^{\varrho}) = \int_{\mathbb{T}} \mu(d\theta) \left( \frac{|b_{\mu,\alpha}(\theta)|^{1/(1-\varrho)}}{g(\theta)} \right)^{1-\varrho} g(\theta) \geq \left[ \mu\left(|b_{\mu,\alpha}|^{1/(1-\varrho)}\right) \right]^{1-\varrho} \quad (90) $$

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_{\mu,\alpha}|^{1/(1-\varrho)}$. ■

The next lemma shows that this choice leads to a non-positive second term, thus implying that $\Psi_{\alpha}(\zeta k) \leq \Psi_{\alpha}(\mu k)$.

**Lemma 5** Assume (A1). Let $\alpha \in \mathbb{R} \setminus \{1\}$ and assume that $\varrho$ satisfies (A2). Then $\eta = (1 - \alpha)/(1 - \varrho)$ satisfies (88) for any $\mu \in M_1(\mathbb{T})$ such that $\mu(|b_{\mu,\alpha}|) < \infty$.

**Proof** We apply (90) with $g = 1$ and get that

$$ [\mu (|b_{\mu,\alpha}|^{1/(1-\varrho)})]^{1-\varrho} \leq \mu (|b_{\mu,\alpha}|) < \infty. \quad (91) $$

Then (88) can be readily checked with $\eta = (1 - \alpha)/(1 - \varrho)$. Set $\phi = \eta/(1 - \alpha)$. Using that $\mu(|b_{\mu,\alpha}|) < \infty$ when $\phi < 0$ and (A1) for $\phi > 0$, we obtain $\mu(|b_{\mu,\alpha}|^{\phi}) > 0$, which concludes the proof. ■

While Lemma 5 seems to advocate for $g = d\zeta/d\mu$ to be proportional to $|b_{\mu,\alpha}|^{1/(1-\varrho)}$, notice that this choice of $g$ might not be optimal to minimize $\Psi_{\alpha}(\zeta k) - \Psi_{\alpha}(\mu k)$, as $A(\mu, \zeta)$ also depends on $g$ through $\zeta$. In the next lemma, we thus propose another choice of the tuning parameter $\eta$, which also satisfies (88) for any $\mu \in M_1(\mathbb{T})$ such that $\mu(|b_{\mu,\alpha}|) < \infty$.

**Lemma 6** Assume (A1). Let $\alpha \in \mathbb{R} \setminus \{1\}$ and assume that $\varrho$ satisfies (A2). Let $\mu \in M_1(\mathbb{T})$ be such that $\mu(|b_{\mu,\alpha}|) < \infty$. Further assume that $|\varrho| \geq 1$, then $\eta = (\alpha - 1)/\varrho$ satisfies (88).
**Proof** Setting $g \propto |b_{\mu,\alpha}|^{-1/\rho}$, we get

$$ \mu(|b_{\mu,\alpha}|g^{\rho}) = \mu(|b_{\mu,\alpha}|^{1-\rho/\rho})[\mu(|b_{\mu,\alpha}|^{-1/\rho})]^{-\rho} = [\mu(|b_{\mu,\alpha}|^{-1/\rho})]^{-\rho} \le \mu(|b_{\mu,\alpha}|) $$

where the last inequality follows from Jensen's inequality applied to the convex function $u \mapsto u^{-\rho}$ (since $|\rho| \ge 1$). Since $\mu(|b_{\mu,\alpha}|) < \infty$, the parameter $\eta = (\alpha - 1)/\rho$ satisfies (88). Set $\phi = \eta/(1 - \alpha)$. Using that $\mu(|b_{\mu,\alpha}|) < \infty$ when $\phi < 0$ and (A1) for $\phi > 0$, we obtain $\mu(|b_{\mu,\alpha}|^{\phi}) > 0$, which concludes the proof. ■

Lemma 5 and Lemma 6 allow us to define a range of values for $\eta$ that decreases $\Psi_{\alpha}$ after one transition of the Power Descent, under the assumption that $\rho$ satisfies (A2). Now, in order to prove Proposition 2 and given $\alpha \in \mathbb{R} \setminus \{1\}$, we need to check which values of $\rho$ satisfy the conditions expressed in (A2).

**Proof of Proposition 2** The proof consists in verifying that we can apply Proposition 6, that is, given $\alpha \in \mathbb{R} \setminus \{1\}$, we must find a range of constants $\rho$ which satisfy (A2). We then use Lemma 5 or Lemma 6 to deduce that, for the provided constants $\eta$, (88) holds.

(i) Assumption (A2) holds for all $\rho < 0$, with $f_{\alpha,\rho}(u) = -\log(u)/\rho$. Moreover, by definition of $b_{\mu,\alpha}$, we get for all $n \ge 1$,

$$ \mu(|b_{\mu,\alpha}|) = \int_{\mathbb{Y}} \mu k(y) \frac{p(y)}{\mu k(y)} \nu(dy) = \int_{\mathbb{Y}} p(y) \nu(dy) < \infty. $$

Combining with Lemma 5 and Lemma 6, (88) holds for all $\mu \in M_1(\mathbb{T})$ and for any $\eta \in (0, 1]$.

(ii) Observing that for $\alpha \notin \{0, 1\}$,

$$ f_{\alpha,\rho}(u) = \frac{1}{\alpha(\alpha - 1)} (u^{\alpha/\rho} - 1), $$

we get that (A2) holds for $\rho \le \alpha$ if $\alpha < 0$ Lemmas 5 and 6 provide the corresponding ranges for $\eta$ in Cases (i) and (ii). To finish the proof, we now show that for all $\mu \in M_1(\mathbb{T})$, $\mu(|b_{\mu,\alpha}|)$ is finite, so that Lemmas 5 and 6 can indeed be applied.

Since $uf'_{\alpha}(u) = \alpha f_{\alpha}(u) + 1/(\alpha - 1)$, we have, for all $n \ge 1$,

$$ \begin{aligned} \mu(|b_{\mu,\alpha}|) &= \int_{\mathbb{Y}} \left| \left( \frac{\mu k(y)}{p(y)} \right) f'_{\alpha} \left( \frac{\mu k(y)}{p(y)} \right) \right| p(y) \nu(dy) && (92) \\ &\le |\alpha| \int_{\mathbb{Y}} \left| f_{\alpha} \left( \frac{\mu k(y)}{p(y)} \right) \right| p(y) \nu(dy) + \frac{1}{|\alpha - 1|} \end{aligned} $$

Using that $\Psi_{\alpha}(\mu k) > -\infty$ (which is a consequence of (A1) and of Jensen's inequality applied to the convex function $u \mapsto uf_{\alpha}(1/u)$), the r.h.s is finite if and only if $\Psi_{\alpha}(\mu k)$ is finite, which is what we have assumed and thus the proof is finished. ■
## C. Deferred Proofs and Results of Section 6

We start with a useful lemma.

### C.1 A Useful Lemma

The following lemma will be used for the proofs of Theorems 4 and 5 and Proposition 3.

**Lemma 7** Let $k^{(o)}$ satisfy (B1). We have, for all $\zeta, \zeta' \in \text{Int} E_0$,

$$ \log \left( \frac{k^{(o)}(\zeta', y)}{k^{(o)}(\zeta, y)} \right) \geq \langle S(y) - \nabla A(\zeta'), \zeta' - \zeta \rangle_E, \quad \text{with equality if and only if } \zeta' = \zeta, \quad (93) $$

$$ 0 \geq \int_{\mathcal{Y}} \log \left( \frac{k^{(o)}(\zeta', y)}{k^{(o)}(\zeta, y)} \right) k^{(o)}(\zeta, y) \nu(dy) \geq \langle \nabla A(\zeta) - \nabla A(\zeta'), \zeta' - \zeta \rangle_E. \quad (94) $$

**Proof** Let $\zeta, \zeta' \in \text{Int} E_0$. We have, since $\text{Int} E_0$ is convex,

$$ \begin{aligned} \log \left( \frac{k^{(o)}(\zeta', y)}{k^{(o)}(\zeta, y)} \right) &= \langle S(y), \zeta' - \zeta \rangle_E + A(\zeta) - A(\zeta') \\ &= \langle S(y), \zeta' - \zeta \rangle_E + \langle \nabla A(\zeta'), \zeta - \zeta' \rangle_E \\ &\quad + (\zeta - \zeta')^T \left( \int_0^1 (1-t) \left[ \nabla \nabla^T A(t\zeta + (1-t)\zeta') \right] dt \right) (\zeta - \zeta'). \end{aligned} $$

The inequality and the equality case of (93) then follow from the fact that the hessian $\nabla \nabla^T A$ is positive definite in $\text{Int} E_0$. Using (93) and (32), we further obtain the lower bound displayed in (94) on the negated exclusive Kullback-Leibler divergence between the two (absolutely continuous w.r.t. $\nu$) probability distributions with probability density functions $k^{(o)}(\zeta, \cdot)$ and $k^{(o)}(\zeta', \cdot)$ respectively. ■

### C.2 Proofs of Section 6.2

The following lemma will be useful.

**Lemma 8** Let $\text{let } p$ and $g$ be functions measurable from $(\mathcal{Y}, \mathcal{Y})$ to the Borel sets of $\mathbb{R}_+$ such that

$$ 0 < \int_{\mathcal{Y}} p(y) (1 + g(y))^{1/(1-\alpha)} \nu(dy) < \infty. \quad (95) $$

Let $\alpha \in [0, 1)$ and $k: \mathbb{T} \times \mathcal{Y} \to \mathbb{R}_+$ be a positive kernel density with respect to $\nu$. Let $\theta_0 \in \mathbb{T}$ and $\mu$ be a probability on $\mathbb{T}$ such that

$$ \inf_{y \in \mathcal{Y}} \frac{k(\theta_0, y)}{\mu k(y)} > 0. \quad (96) $$

Then, setting

$$ \varphi(y) = k(\theta_0, y)^\alpha \left( \frac{p(y)}{\mu k(y)} \right)^{1-\alpha}, $$
we have $\int \varphi \, d\nu > 0$ and the probability density function $\tilde{\varphi} = \varphi / (\int \varphi \, d\nu)$ satisfies

$$ \int \tilde{\varphi}(y) g(y) \nu(dy) < \infty. \qquad (97) $$

**Proof** We have $\mu k > 0$ by assumption on $k$ and $\mu$. Thus $\varphi$ is well defined and $\int \varphi \, d\nu > 0$ as a consequence of $p \ge 0$ with the left-hand side of (95). We now check (97) in which $\tilde{\varphi}$ can equivalently be replaced by $\varphi$. We have

$$ \begin{aligned} \int \varphi(y) g(y) \nu(dy) &= \int k(\theta_0, y) \left( \frac{p(y)}{\mu k(y)} \right)^{1-\alpha} g(y) \nu(dy) \\ &\le \left( \int k(\theta_0, y) \frac{p(y)}{\mu k(y)} g(y)^{1/(1-\alpha)} \nu(dy) \right)^{\alpha}, \end{aligned} $$

by Jensen's inequality and concavity of $x \mapsto x^{1-\alpha}$ for $\alpha \in [0, 1)$. Hence we obtain (97) as a consequence of (95) and (96). ■

**Proof of Theorem 4** Let $\zeta_0 \in \mathcal{O}$ and $b \ge 0$. In the following, we denote by $C_0(\zeta)$ the value of the integral in the argmax in (35). Using that $\tilde{\varphi}$ is a probability density function paired up with (32) and (34), we have

$$ \begin{aligned} \int_{\mathcal{Y}} [\tilde{\varphi}(y) + b k^{(o)}(\zeta_0, y)] \nu(dy) &= 1 + b > 0, \\ \int_{\mathcal{Y}} [\tilde{\varphi}(y) + b k^{(o)}(\zeta_0, y)] S(y) \nu(dy) &= \int_{\mathcal{Y}} \tilde{\varphi} S \, d\nu + b \nabla A(\zeta_0) \in E. \end{aligned} $$

By (93) in Lemma 7, for all $\zeta' \ne \zeta$ in $\mathcal{O}$ and all $y \in \mathcal{Y}$, we have

$$ \log \left( \frac{k^{(o)}(\zeta', y)}{k^{(o)}(\zeta, y)} \right) > \langle S(y) - \nabla A(\zeta'), \zeta' - \zeta \rangle_E . $$

Hence, combining with the three previous identities, we get that, for all $\zeta' \ne \zeta$ in $\mathcal{O}$,

$$ \begin{aligned} C_0(\zeta') - C_0(\zeta) &= \int_{\mathcal{Y}} [\tilde{\varphi}(y) + b k^{(o)}(\zeta_0, y)] \log \left( \frac{k^{(o)}(\zeta', y)}{k^{(o)}(\zeta, y)} \right) \nu(dy) \\ &> \left\langle \int_{\mathcal{Y}} \tilde{\varphi} S \, d\nu + b \nabla A(\zeta_0) - (1+b) \nabla A(\zeta'), \zeta' - \zeta \right\rangle_E \\ &= (1+b) \langle \mathbf{s}^* - \nabla A(\zeta'), \zeta' - \zeta \rangle_E, \end{aligned} \qquad (98) $$

where we used the definition of $\mathbf{s}^*$ in (36). Now, since $\nabla A$ is a $C^\infty$-diffeomorphism on $\text{Int} E_0$, if $\mathbf{s}^*$ belongs to $\nabla A(\mathcal{O})$ then there exists a unique $\zeta^* \in \mathcal{O}$ such that $\nabla A(\zeta^*) = \mathbf{s}^*$. As a result, plugging this into (98) and setting $\zeta' = \zeta^*$, we have that $C_0(\zeta^*) > C_0(\zeta)$ for all $\zeta \ne \zeta^*$. This shows that $\zeta^*$ is the (unique) solution to the argmax problem (35).

We conclude this proof with the proof of the reciprocal implication. Suppose that the argmax problem (35) has a solution $\zeta \in \mathcal{O}$ so that, for all $\zeta' \in \mathcal{O}$, $C_0(\zeta') - C_0(\zeta) \le 0$. Using (98), this would imply $\langle \mathbf{s}^* - \nabla A(\zeta'), \zeta' - \zeta \rangle < 0$ for all $\zeta' \in \mathcal{O} \setminus \{\zeta\}$. Since $\mathcal{O}$ is an open set, we can take $\zeta' = \zeta + \epsilon (\mathbf{s}^* - \nabla A(\zeta'))$ with $\epsilon > 0$ small enough, which gives

$$ \langle \mathbf{s}^* - \nabla A(\zeta'), \mathbf{s}^* - \nabla A(\zeta) \rangle < 0. $$
Now letting $\epsilon \downarrow 0$ and since $\nabla A$ is continuous, we get that $s^* - \nabla A(\zeta) = 0$ and thus that $s^*$ belongs to $\nabla A(\mathcal{O})$. ■

**Proof of Corollary 6** First note that letting $\nu$ be the $d$-dimensional Lebesgue measure and setting $S(y) = (y, -yy^T/2)$, the condition (34) written in Theorem 4 takes the form

$$ \int_{\mathbb{Y}} \|y\|^2 \check{\varphi}(y) \, dy < \infty. \qquad (99) $$

Furthermore, setting $v(\theta) = (\Sigma^{-1}m, \Sigma^{-1})$ with $\theta = (m, \Sigma)$, we obtain that $v(\mathbb{T}) = \mathbb{R}^d \times \mathcal{M}_{>0}(d)$, which is the whole interior set $\text{Int}E_0$, seen as an open set of $E = \mathbb{R}^d \times \mathbb{R}^{d \times d}$. Hence Theorem 4 applies and we can use the interpretation (38) of the equation (37) with $k(\theta^*, y) = \mathcal{N}(y; m^*, \Sigma^*)$ and $k(\theta_0, y) = \mathcal{N}(y; m_0, \Sigma_0)$, that is

$$ \int_{\mathbb{Y}} S(y) \mathcal{N}(y; m^*, \Sigma^*) \nu(dy) = \int_{\mathbb{Y}} S(y) \psi(y) \nu(dy), \qquad (100) $$

where $\psi$ is the mixture

$$ \psi(y) = \frac{1}{1+b} \check{\varphi}(y) + \frac{b}{1+b} \mathcal{N}(y; m_0, \Sigma_0), \quad y \in \mathbb{Y}. $$

Using now that $S(y) = (y, -yy^T/2)$, (100) is equivalent to saying that the distributions on both sides of this equation have the same mean and covariance matrix. Since the left-hand side distribution is a Gaussian distribution, we directly obtain that the solution of (100) is given by $\theta^* = (m^*, \Sigma^*)$ with $m^*$ and $\Sigma^*$ being the mean and covariance matrix of the random variable $Y$ with density $\psi$.

Observe that these mean and covariance matrix are well-defined due to the condition (99). This means that we only have to check that the covariance matrix is positive definite in order to have a solution in $\mathbb{T} = \mathbb{R}^d \times \mathcal{M}_{>0}(d)$. If $b > 0$, then $1/(1+b) < 1$ and the fact that $\theta \in \mathbb{T}$ guarantees that the covariance matrix of $Y$ is positive definite. If $b = 0$, then $Y$ has density $\check{\varphi}$ and the non-degenerate condition on $\check{\varphi}$ guarantees that the covariance is also positive definite in this case. This concludes the proof. ■

**Proof of Corollary 7** Recall from (16) that $\varphi_{j,n}^{(\alpha)}(y) = k(\theta_{j,n}, y)(p(y)/\mu_n k(y))^{1-\alpha}$ with $k(\theta_{j,n}, y) > 0$. Then Lemma 8 gives that

$$ 0 < \int_{\mathbb{Y}} (1 + \|S(y)\|_E) \varphi_{j,n}^{(\alpha)}(y) \nu(dy) < \infty. \qquad (101) $$

We can thus apply Theorem 4 with $\zeta_0 = v(\theta_{j,n})$, $\mathcal{O} = v(\mathbb{T})$, $b = b_{j,n} / \int_{\mathbb{Y}} \varphi_{j,n}^{(\alpha)} d\nu$ and $\check{\varphi} = \varphi_{j,n}^{(\alpha)} / \int \varphi_{j,n}^{(\alpha)} d\nu$, which satisfies (34). ■

### C.3 Proofs of Section 6.3

**Proof of Proposition 3** Following the proof of Theorem 4 but with $b = 0$, we have using (98) that: for all $\zeta', \zeta \in \text{Int}E_0$,

$$ \left\langle \nabla A(\zeta) - \int_{\mathbb{Y}} S \check{\varphi} \, d\nu, \zeta' - \zeta \right\rangle_E \le g^{(o)}(\zeta') - g^{(o)}(\zeta) \le \left\langle \nabla A(\zeta') - \int_{\mathbb{Y}} S \check{\varphi} \, d\nu, \zeta' - \zeta \right\rangle_E. $$
Since $\nabla A$ is continuous in $\text{Int} E_0$, we get that $g^{(0)}$ is $C^1$ over $\text{Int} E_0$ and, for all $\zeta \in \text{Int} E_0$,

$$ \nabla g^{(0)}(\zeta) = \nabla A(\zeta) - \int_{\mathbb{Y}} S \varphi \, d\nu . $$

The $\beta_0$-smoothness condition follows from the mean value theorem applied to $\nabla g^{(0)}$ whose Jacobian is $\nabla \nabla^T A$. ■

**Proof of Corollary 8** The context is similar to the proof of Corollary 6 except that the covariance matrix is known and thus the sufficient statistic reduces to $S(y) = y$ and the canonical parameter to $v(\theta) = \Sigma^{-1}\theta$. The given expressions of $h$ and $A$ follow. Then we have, for all $\zeta \in \mathbb{R}^d$, $\nabla A(\zeta) = \Sigma\zeta$. Applying Proposition 3, we get that

$$ \nabla g^{(0)}(\zeta) = \Sigma\zeta - \int_{\mathbb{Y}} y\varphi(y) \, dy . $$

This gives that $g^{(0)}$ is $\beta_0$-smooth over $\mathbb{R}^d$ with $\beta_0$ equal to the maximal eigenvalue of $\Sigma$ and leads to the gradient step

$$ \zeta = \zeta_0 - \frac{\gamma}{\beta_0} \left( \Sigma\zeta_0 - \int_{\mathbb{Y}} y\varphi(y) \, dy \right) . $$

Multiplying on both sides by $\Sigma$ leads to the gradient-based update (46).

We now derive (47). Using that $g = g^{(0)} \circ v$, we get that $\nabla g(\theta) = \Sigma^{-1} \nabla g^{(0)}(\Sigma^{-1}\theta)$ and thus

$$ \nabla g(\theta) = \Sigma^{-1} \left( \theta - \int_{\mathbb{Y}} y \, \check{\varphi}(y) \right) . $$

We obtain that $g$ is $\beta$-smooth over $\mathbb{R}^d$ with $\beta$ equal to the maximal eigenvalue of $\Sigma^{-1}$, which leads to the gradient update (47). ■

## C.4 Proofs of Section 6.4

### C.4.1 PROOFS OF PRELIMINARY RESULTS

It will be convenient to use the shorthand notation: for all $(\ell, y, \xi) \in \mathbb{Z} \times \mathbb{Y} \times F$,

$$ k_{\xi}(\ell, y) := k^{(0)}(\ell(\xi), y), $$

$$ \check{k}_{\xi, \tau}(\ell, y) := \frac{k_{\xi}(\ell, y)}{\tau k_{\xi}(y)}, $$

where $\tau k_{\xi}$ denotes the mapping $y \mapsto \int k_{\xi}(\ell, y) \tau(d\ell)$.
**Proof of Proposition 4** We have, for all $\tau_0, \tau \in M, \xi_0, \xi \in F$ and $y \in Y$, using the definitions above and that of $D_{\tau, \tau_0}$,

$$ 
\begin{align*}
\log \frac{\tau k_{\xi}(y)}{\tau_0 k_{\xi_0}(y)} &= \log \int_Z \left( \frac{k_{\xi}(\ell, y)}{k_{\xi_0}(\ell, y)} \right) \check{k}_{\xi_0, \tau_0}(\ell, y) \tau(d\ell) \\
&= \log \int_Z \left( \frac{k_{\xi}(\ell, y)}{k_{\xi_0}(\ell, y)} D_{\tau, \tau_0}(\ell) \right) \check{k}_{\xi_0, \tau_0}(\ell, y) \tau_0(d\ell) \\
&\geq \int_Z \check{k}_{\xi_0, \tau_0}(\ell, y) \log \left( \frac{k_{\xi}(\ell, y)}{k_{\xi_0}(\ell, y)} D_{\tau, \tau_0}(\ell) \right) \tau_0(d\ell),
\end{align*}
$$

where we used the concavity of the log function and the fact that $\check{k}_{\xi_0, \tau_0}(\cdot, y)$ is a probability density function with respect to $\tau_0$. Let $\xi = \Upsilon(\vartheta)$ and $\xi_0 = \Upsilon(\vartheta_0)$ and set $\theta = (\vartheta, \tau)$ and $\theta_0 = (\vartheta_0, \tau_0)$, so that, by Definition 2, we have $\tau k_{\xi}(y) = k^{(1)}(\theta, y)$ and $\tau_0 k_{\xi_0}(y) = k^{(1)}(\theta_0, y)$. Applying the previous bound and integrating w.r.t. $\check{\varphi}(y)\nu(dy)$ thus yields

$$ 
\begin{align*}
\int_Y \check{\varphi}(y) \log \left( \frac{k^{(1)}(\theta, y)}{k^{(1)}(\theta_0, y)} \right) \nu(dy) \\
&\geq \int_Y \left( \int_Z \varphi(\ell, y) \log \left( \frac{k_{\xi}(\ell, y)}{k_{\xi_0}(\ell, y)} D_{\tau, \tau_0}(\ell) \right) \tau_0(d\ell) \right) \nu(dy), \quad (102)
\end{align*}
$$

where we used that $\check{\varphi}(y)\check{k}_{\xi_0, \tau_0}(\ell, y) = \varphi(\ell, y)$ as a consequence of (52). We conclude by noting that (50), (51) and (102) imply (49). $\blacksquare$

**Proof of Corollary 9** Suppose that $\Upsilon(\vartheta)$ belongs to the argmax (54), and let us show that (51) holds. Note that the objective function in the argmax is zero for $\xi = \xi_0$ and thus must be non-negative for $\xi = \Upsilon(\vartheta)$ in the argmax. Note also that the left-hand side of (51) is the sum of this objective function at $\xi = \Upsilon(\vartheta)$ with the term

$$ - \int_Y \left( \int_Z b_0(\ell) k_{\xi_0}(\ell, y) \log \left( \frac{k_{\Upsilon(\vartheta)}(\ell, y)}{k_{\xi_0}(\ell, y)} \right) \tau_0(d\ell) \right) \nu(dy), \quad (103) $$

and it only remains to show that this term is also non-negative. To conclude, we observe that, for $\tau_0$-a.e. $\ell \in Z$, Condition (ii) of Definition 1 implies that $k_{\xi_0}(\ell, \cdot)$ and $k_{\Upsilon(\vartheta)}(\ell, \cdot)$ are probability densities, so that

$$ - \int_Y k_{\xi_0}(\ell, y) \log \left( \frac{k_{\Upsilon(\vartheta)}(\ell, y)}{k_{\xi_0}(\ell, y)} \right) \nu(dy) \geq 0, $$

and thus, $b_0$ being non-negative,

$$ - \int_Z \left( \int_Y b_0(\ell) k_{\xi_0}(\ell, y) \log \left( \frac{k_{\Upsilon(\vartheta)}(\ell, y)}{k_{\xi_0}(\ell, y)} \right) \nu(dy) \right) \tau_0(d\ell) \geq 0. $$

Using that $\log_+(x) \leq x$ for all $x \geq 0$, the positive part inside this double integral is at most $b_0(\ell) k_{\Upsilon(\vartheta)}(\ell, y)$, whose integral w.r.t. $\nu(dy)\tau_0(d\ell)$ is $\int b_0 d\tau_0 < \infty$ by assumption. Therefore,
the same inequality holds when reversing the order of integration in the previous display. We get that the expression in (103) is non-negative, which concludes the proof. ■

The following lemma is used to illustrate (C1).

**Lemma 9** Consider an ELM family satisfying (C1). Then, for all $\vartheta \in \tilde{T}$, $\tau \in M$ and $\tau$-a.e. $\ell$, $\ell \circ \Upsilon(\tilde{T})$ is an open subset included in $\text{Int}E_0$.

**Proof** In (C1), the linear operator $\ell \in Z$ being full rank implies that the image of an open set is an open subset of $E$. Indeed, take $F' \subset F$ a linear supplementary space of the kernel of $\ell$. Then $\ell|_{F'}$ is a bijective linear mapping from $F'$ to $E$ and thus, for any open set $O$ included in $F$, $\ell(O) = \ell(O \cap F') = \ell|_{F'}(O \cap F')$ is an open subset of $E$ since $O \cap F'$ is an open subset of $F'$ and $\ell|_{F'}$ is bicontinuous. Hence (C1) implies that $\ell \circ \Upsilon(\tilde{T})$ is an open subset of $E$ for all $\tau \in M$ and $\tau$-almost all $\ell \in Z$, which, by Assumption (ii) from Definition 2, must be included in $\text{Int}E_0$. ■

We conclude with the following lemma, which relates Assumption (C4) to a kind of identifiability property.

**Lemma 10** Consider an ELM family and let $\mathbb{P}_{\vartheta, \tau}$ denote the probability of $(Y, L)$ defined on $Y \times Z$ as above. Then (C4) is equivalent to have that for all $(\vartheta, \tau)$, $(\vartheta', \tau') \in \tilde{T} \times M$, $\mathbb{P}_{\vartheta, \tau} = \mathbb{P}_{\vartheta', \tau'}$ if and only if $\Upsilon(\vartheta) = \Upsilon(\vartheta')$ and $\tau = \tau'$.

**Proof** Let $(\vartheta, \tau)$, $(\vartheta', \tau') \in \tilde{T} \times M$. Note that $\Upsilon(\vartheta) = \Upsilon(\vartheta')$ and $\tau = \tau'$ is always sufficient to have $\mathbb{P}_{\vartheta, \tau} = \mathbb{P}_{\vartheta', \tau'}$. Also, since $\tau$ is the marginal over $Z$ of $\mathbb{P}_{\vartheta, \tau}$ and $k^{(o)}(\ell \circ \Upsilon(\vartheta), \cdot)$ is the conditional density of $Y$ given $L = \ell$ with respect to $\nu$, it is clear that $\mathbb{P}_{\vartheta, \tau} = \mathbb{P}_{\vartheta', \tau'}$ is equivalent to have $\tau = \tau'$ and

$$k^{(o)}(\ell \circ \Upsilon(\vartheta), y) = k^{(o)}(\ell \circ \Upsilon(\vartheta'), y) \quad \text{for } \nu\text{-a.e. } y \text{ and } \tau\text{-a.e. } \ell. \quad (104)$$

By (B1), which is assumed for an ELM family in Definition 2, for all $\zeta, \zeta' \in E_0$, we have $k^{(o)}(\zeta, y) = k^{(o)}(\zeta', y)$ for $\nu$-a.e. $y$ if and only if $\zeta = \zeta'$. By (ii) in Definition 2, we have $\ell \circ \Upsilon(\vartheta) \in E_0$ for all $\vartheta \in \tilde{T}$, $\tau \in M$ and $\tau$-a.e. $\ell$. Thus, for all $(\vartheta, \vartheta', \tau) \in \tilde{T}^2 \times M$, Eq. (104) is equivalent to have

$$\ell \circ \Upsilon(\vartheta) = \ell \circ \Upsilon(\vartheta') \quad \text{for } \tau\text{-a.e. } \ell. \quad (105)$$

Hence, to prove the result, we only need to show that (C4) is equivalent to have that, for all $(\vartheta, \vartheta', \tau) \in \tilde{T}^2 \times M$, Eq. (105) implies $\Upsilon(\vartheta) = \Upsilon(\vartheta')$. In fact (C4) exactly says that for all $\vartheta, \vartheta' \in \tilde{T}$, $\Upsilon(\vartheta) \neq \Upsilon(\vartheta')$ implies the opposite of (105). ■

### C.4.2 PROOF OF THEOREM 5

The proof of Theorem 5 relies on a general proposition (Proposition 7) and two technical lemmas (Lemmas 11 and 12). The following definition will be useful to state the proposition.

**Definition 4 (Directional continuity)** Let $g$ be defined on an open subset $O$ of the Euclidean space $F$ and valued in a topological space. Let $\xi \in O$ and $x \in F$ such that $\|x\|_F = 1$. We say that $g$ is continuous at $\xi$ in direction $x$ if $\lim_{t \downarrow 0} g(\xi + tx) = g(\xi)$. ■
**Proposition 7** Consider an ELM family as in Definition 2 and $\Upsilon$ satisfying (C1). Let $\tilde{\varphi}$ be a probability density function w.r.t. $\nu$. Let $\tau_0 \in M$ and $\varphi: \mathbb{Z} \times \mathbb{Y} \to \mathbb{R}_+$ be measurable function such that

$$ \int_{\mathbb{Y} \times \mathbb{Z}} \varphi(\ell, y) \|\ell\|_{\text{op}} \|S(y)\|_{E} \nu(dy) \tau_0(d\ell) < \infty \quad (106) $$

$$ \int_{\mathbb{Z}} \tilde{\varphi}(\ell) \left\| \ell^{\nabla}(\xi) \right\|_{F} \tau_{0}(d\ell) < \infty \quad \text{for all } \xi \in \Upsilon(\tilde{\mathcal{T}}), \quad (107) $$

where we set $\varphi(\ell, y)$ and $\tilde{\varphi}(\ell)$ as in (52) and (59). Define $w_{\tau_0}^{\nabla}: \Upsilon(\tilde{\mathcal{T}}) \to F$ as in (60). Let $b \in \mathbb{R}_+$ and $\xi_0$ in $\Upsilon(\tilde{\mathcal{T}})$ and set $\mathbf{w}^*$ as in (63). Then the reciprocal set

$$ \Xi^* = \{ \xi \in \Upsilon(\tilde{\mathcal{T}}) : w_{\tau_0}^{\nabla}(\xi) = \mathbf{w}^* \} $$

is included in the set of solutions of the argmax problem

$$ \underset{\xi \in \Upsilon(\tilde{\mathcal{T}})}{\operatorname{argmax}} \int_{\mathbb{Y}} \left( \int_{\mathbb{Z}} \left[ \varphi(\ell, y) + b \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \right] \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \tau_0(d\ell) \right) \nu(dy). \quad (108) $$

Moreover the two following assertions hold.

(i) If $w_{\tau_0}^{\nabla}$ is continuous on $\Upsilon(\tilde{\mathcal{T}})$ in every direction, the opposite inclusion is true, namely, the solution set of the argmax problem (108) is included in $\Xi^*$.

(ii) If we have, for all $\xi \neq \xi' \in \Upsilon(\tilde{\mathcal{T}})$,

$$ \int_{\{\ell \in \mathbb{Z} : \ell(\xi - \xi') \neq 0\}} \tilde{\varphi}(\ell) \tau_0(d\ell) > 0, \quad (109) $$

then $\xi \mapsto w_{\tau_0}^{\nabla}(\xi)$ is one-to-one on $\Upsilon(\tilde{\mathcal{T}})$ and thus $\Xi^*$ is either a singleton or empty.

**Proof** In this proof, we take $b \ge 0$ and $\xi_0 \in \Upsilon(\tilde{\mathcal{T}})$, and we denote the integral defining the objective function in the argmax problem (108) by $C(\xi)$.

A first remark is that, by (106), $\int \varphi(\ell, y) \ell^T \circ S(y) \nu(dy) \tau_0(d\ell)$ is well-defined in $F$ and $w_{\tau_0}^{\nabla}(\xi)$ in (60) is well-defined in $F$ for all $\xi \in \Upsilon(\tilde{\mathcal{T}})$ as a consequence of (107). In particular we get that $\mathbf{w}^*$ in (63) is well-defined in $F$ (but does not necessarily belong to $w_{\tau_0}^{\nabla} \circ \Upsilon(\tilde{\mathcal{T}})$).

Next we show that the objective function $C(\xi)$ is well defined and finite for all $\xi \in \Upsilon(\tilde{\mathcal{T}})$ and that Fubini's theorem applies (so that the double integral in (108) can be computed in the opposite order or as a single integral over $\mathbb{Y} \times \mathbb{Z}$ w.r.t. the product measure $\nu \otimes \tau_0$). To this end, we show successively that

$$ \int_{\mathbb{Y} \times \mathbb{Z}} \varphi(\ell, y) \left| \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \right| \nu(dy) \tau_0(d\ell) < \infty, \quad (110) $$

$$ \int_{\mathbb{Y} \times \mathbb{Z}} \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \log_+(\frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)}) \nu(dy) \tau_0(d\ell) < \infty, \quad (111) $$

$$ -\infty < \int_{\mathbb{Z}} \tilde{\varphi}(\ell) \left( \int_{\mathbb{Y}} k^{(o)}(\ell(\xi_0), y) \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \nu(dy) \right) \tau_0(d\ell) \le 0, \quad (112) $$

$$ \int_{\mathbb{Y} \times \mathbb{Z}} \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \left| \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \right| \nu(dy) \tau_0(d\ell) < \infty. \quad (113) $$
Let us now show that (110) holds. Let $\ell \in \mathbb{Z}$ such that $\ell \circ \Upsilon(\tilde{T}) \subseteq \text{Int}E_0$, which by Condition (ii) and (C1), is true for $\tau_0$-a.e. $\ell \in \mathbb{Z}$. Applying (93) in Lemma 7 with $(\zeta, \zeta') = (\ell(\xi_0), \ell(\xi))$ as well as $(\zeta, \zeta') = (\ell(\xi), \ell(\xi_0))$ and using the definition of $\ell^T$ in (55) and that of $\ell^\nabla$ in (56), we get that

$$ \langle \ell^T \circ S(y) - \ell^\nabla(\xi), \xi - \xi_0 \rangle_F \le \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \le \langle \ell^T \circ S(y) - \ell^\nabla(\xi_0), \xi - \xi_0 \rangle_F. $$

It follows that, for all $\xi_0, \xi \in \Upsilon(\tilde{T})$,

$$ \left| \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \right| \le \left( \left\| \ell^T \circ S(y) \right\|_F + \max \left( \left\| \ell^\nabla(\xi_0) \right\|_F, \left\| \ell^\nabla(\xi) \right\|_F \right) \right) \|\xi - \xi_0\|_F . \quad (114) $$

Using $\|\ell^T \circ S(y)\|_F \le \|\ell^T\|_{\text{op}} \|S(y)\|_E = \|\ell\|_{\text{op}} \|S(y)\|_E$, (106) and (107), we obtain (110).

To show (111), we use the bound $\log_+(x) \le x$ for all $x \ge 0$ so that the left-hand side of (111) is upper-bounded by $\int_{\mathbb{Y} \times \mathbb{Z}} \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi), y) \nu(dy)\tau_0(d\ell) = \int_{\mathbb{Z}} \tilde{\varphi}(\ell) \tau_0(d\ell) = 1$.

Next, to show (112), we observe that, as a consequence of (94), (55) and (56), we have

$$ \begin{aligned} 0 \ge \int_{\mathbb{Y}} k^{(o)}(\ell(\xi), y) \log \left( \frac{k^{(o)}(\ell(\xi), y)}{k^{(o)}(\ell(\xi_0), y)} \right) \nu(dy) &\ge \langle \ell^\nabla(\xi_0) - \ell^\nabla(\xi), \xi - \xi_0 \rangle_F \\ &\ge - \| \ell^\nabla(\xi_0) - \ell^\nabla(\xi) \|_F \|\xi - \xi_0\|_F . \end{aligned} $$

Hence, using (107), we get (112).

Finally, Tonelli's theorem gives us that the same positive part as in (111) has a finite integral when integrating in the same order as in (112), so that (111) and (112) imply (113).

We now prove the inclusion of $\Xi^*$ in the set of argmax solutions. Since Fubini's theorem applies in the definition of the objective function, we can write, for all $\xi, \xi' \in \Upsilon(\tilde{T})$,

$$ C(\xi') - C(\xi) = \int_{\mathbb{Y} \times \mathbb{Z}} \left[ \varphi(\ell, y) + b \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \right] \log \left( \frac{k^{(o)}(\ell(\xi'), y)}{k^{(o)}(\ell(\xi), y)} \right) \nu(dy)\tau_0(d\ell). $$

Using (93) in Lemma 7 this time with $\zeta = \ell(\xi)$ and $\zeta' = \ell(\xi')$, we obtain that

$$ C(\xi') - C(\xi) \ge \int_{\mathbb{Y} \times \mathbb{Z}} \left[ \varphi(\ell, y) + b \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \right] \langle \ell^T \circ S(y) - \ell^\nabla(\xi'), \xi' - \xi \rangle_E \nu(dy)\tau_0(d\ell). $$

Using (106), (107) to integrate in any convenient order, we have

$$ \begin{aligned} \int_{\mathbb{Y} \times \mathbb{Z}} \varphi(\ell, y) \langle \ell^T \circ S(y) - \ell^\nabla(\xi'), \xi' - \xi \rangle_E \nu(dy)\tau_0(d\ell) \\ = \left\langle \int_{\mathbb{Y} \times \mathbb{Z}} \varphi(\ell, y) \ell^T \circ S(y) \nu(dy)\tau_0(d\ell) - w_{\tau_0}^\nabla(\xi'), \xi' - \xi \right\rangle_E, \end{aligned} $$

where we have used the definition of $w_{\tau_0}^\nabla$ in (60), and

$$ \begin{aligned} \int_{\mathbb{Y} \times \mathbb{Z}} \tilde{\varphi}(\ell) k^{(o)}(\ell(\xi_0), y) \langle \ell^T \circ S(y) - \ell^\nabla(\xi'), \xi' - \xi \rangle_E \nu(dy)\tau_0(d\ell) \\ = \left\langle w_{\tau_0}^\nabla(\xi_0) - w_{\tau_0}^\nabla(\xi'), \xi' - \xi \right\rangle_E. \end{aligned} $$
where we have used once again the definition of $w_{\tau_0}^{\nabla}$ in (60) as well as the expression of $\nabla A$ given in (32). Plugging these two equalities in the previous inequality, we finally get that, for all $\xi, \xi' \in \Upsilon(\tilde{T})$,

$$ \begin{aligned} C(\xi') - C(\xi) &\geq \left\langle \int_{\Upsilon \times \mathbb{Z}} \varphi(\ell, y) \ell^T \circ S(y) \nu(dy) \tau_0(d\ell) - (1+b) w_{\tau_0}^{\nabla}(\xi') + b w_{\tau_0}^{\nabla}(\xi_0), \xi' - \xi \right\rangle_E \\ &= (1+b) \left\langle \mathbf{w}^* - w_{\tau_0}^{\nabla}(\xi'), \xi' - \xi \right\rangle_E, \end{aligned} \quad (115) $$

Hence any $\xi' \in \Xi^*$ must be a a solution of the argmax (108).

Let us now prove Assertion (i). We use again (115) this time with $\xi$ in the set of the argmax solutions and $\xi'$ in the neighborhood of $\xi$ (remember that $\Upsilon(\tilde{T})$ is assumed to be open) such that $\xi' - \xi = \epsilon (\mathbf{w}^* - w_{\tau_0}^{\nabla}(\xi))$ with $\epsilon > 0$ small enough. We obtain that

$$ \left\langle \mathbf{w}^* - w_{\tau_0}^{\nabla}(\xi'), \mathbf{w}^* - w_{\tau_0}^{\nabla}(\xi) \right\rangle_E \leq 0. $$

The assumption that $w_{\tau_0}^{\nabla}$ is continuous at $\xi$ in direction $(\mathbf{w}^* - w_{\tau_0}^{\nabla}(\xi))$ gives us that $\mathbf{w}^* = w_{\tau_0}^{\nabla}(\xi)$, that is, $\xi \in \Xi^*$ and the opposite inclusion is shown.

We conclude with the proof of Assertion (ii). By (60), for all $\xi, \xi' \in \Upsilon(\tilde{T})$, we can write

$$ w_{\tau_0}^{\nabla}(\xi) - w_{\tau_0}^{\nabla}(\xi') = \int_{\mathbb{Z}} \tilde{\varphi}(\ell) (\ell^{\nabla}(\xi) - \ell^{\nabla}(\xi')) \tau_0(d\ell). $$

Since $\nabla A$ has a positive definite Jacobian on the convex set $\text{Int}E_0$, we have, for all $\zeta \neq \zeta' \in \text{Int}E_0$, $\langle \nabla A(\zeta) - \nabla A(\zeta'), \zeta - \zeta' \rangle_E > 0$.

Letting $\xi \neq \xi' \in \Upsilon(\tilde{T})$ and $\ell \in \mathbb{Z}$ be such that $\ell(\xi), \ell(\xi') \in \text{Int}E_0$ and $\ell(\xi - \xi') \neq 0$, we can apply the previous inequality with $\zeta = \ell(\xi)$ and $\zeta' = \ell(\xi')$ leading to

$$ \left\langle \ell^{\nabla}(\xi) - \ell^{\nabla}(\xi'), \xi - \xi' \right\rangle_F > 0. $$

Combining this result with the condition (109) gives us that for all $\xi \neq \xi' \in \Upsilon(\tilde{T})$

$$ \left\langle w_{\tau_0}^{\nabla}(\xi) - w_{\tau_0}^{\nabla}(\xi'), \xi - \xi' \right\rangle_F > 0. $$

In particular, $w_{\tau_0}^{\nabla}$ is one-to-one and the proof is concluded.

<figure><img src="image_10.png" /></figure>

Let $d_F$ denote the dimension of $F$ and let us now introduce and prove two key lemmas.

**Lemma 11** Let $\xi \in \Upsilon(\tilde{T})$ and let $x \in F$ be such that $\|x\|_F = 1$. Let us denote, for any $\epsilon > 0$,

$$ B_F(\xi, \epsilon) = \{\xi' \in F : \|\xi - \xi'\|_F < \epsilon\}, $$

$$ D_x(\xi, \epsilon) = \{\xi + t \epsilon x : 0 \le t < 1\}. $$

Let $\ell \in \mathbb{Z}$ be such that $\ell \circ \Upsilon(\tilde{T}) \subset \text{Int}E_0$. Then the two following assertions hold.

(i) For any $0 < \epsilon \le \epsilon' < \epsilon''$ such that $D_x(\xi, \epsilon'') \subset \Upsilon(\tilde{T})$, we have

$$ \sup_{\xi' \in D_x(\xi, \epsilon)} |A \circ \ell(\xi') - A \circ \ell(\xi)| \le \epsilon \max \left( \|\ell^{\nabla}(\xi)\|_F, \|\ell^{\nabla}(\xi + \epsilon'x)\|_F \right). \quad (116) $$
(ii) For any $\epsilon > 0$ such that $B_F(\xi, 2(1 + \sqrt{d_F})\epsilon) \subset \Upsilon(\tilde{T})$, we have, for all $\epsilon_0 \in (0, \epsilon]$,

$$ \int_{\Upsilon} e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} \left( \sup_{\xi' \in D_x(\xi, \epsilon_0)} k^{(o)}(\ell(\xi'), y) \right) \nu(dy) \le 2 d_F \max_{\xi' \in \Delta_{x,\epsilon}(\xi)} \left( e^{\epsilon_0' \|\ell^\nabla(\xi')\|_F} \right), \quad (117) $$

where $\epsilon_0' = (1 + 2\sqrt{d_F})\epsilon_0$ and $\Delta_{x,\epsilon}(\xi)$ is a set of $2(d_F + 1)$ points in $B_x(\xi, 2(1 + \sqrt{d_F})\epsilon)$ only depending on $\xi, x$ and $\epsilon$.

**Proof** We prove Assertions (i) and (ii) successively.

**Proof of Assertion (i)** Let $0 < \epsilon \le \epsilon' < \epsilon''$ be such that $D_x(\xi, \epsilon'') \subset \Upsilon(\tilde{T})$. First note that $t \mapsto A \circ \ell(\xi + tx)$ has its first and second derivative on $t \in [0, \epsilon']$ given respectively by

$$ \begin{aligned} h_1 : t &\mapsto \langle \nabla A \circ \ell(\xi + tx), \ell(x) \rangle_E = \langle \ell^\nabla(\xi + tx), x \rangle_F \\ h_2 : t &\mapsto \langle [\nabla \nabla^T A \circ \ell(\xi + tx)] \ell(x), \ell(x) \rangle_E . \end{aligned} $$

Since $\nabla \nabla^T A$ is positive definite on $\mathrm{Int} E_0$, $h_2$ is non-negative, thus $h_1$ is monotonous on $[0, \epsilon']$ and its maximal absolute value is obtained at $t = 0$ or $t = \epsilon'$, meaning that $|h_1(t)| \le \max(|h_1(0)|, |h_1(\epsilon')|)$ for all $t \in [0, \epsilon']$. Therefore, the bound in (116) follows from the mean value theorem combined with the fact that: for all $t \in [0, \epsilon']$ and $\xi' \in D_x(\xi, \epsilon)$,

$$ \begin{aligned} |h_1(t)| &\le \|\ell^\nabla(\xi + tx)\|_F \|x\|_F = \|\ell^\nabla(\xi + tx)\|_F \\ &\le \|\xi' - \xi\|_F \le \epsilon \|x\|_F \le \epsilon \end{aligned} $$

since $\|x\|_F = 1$ by assumption.

**Proof of Assertion (ii)** Take $0 < \epsilon_0 \le \epsilon$ such that $B_x(\xi, 2(1 + \sqrt{d_F})\epsilon) \subset \Upsilon(\tilde{T})$, let $\epsilon_0'$ be as stated and set $\epsilon' = (1 + 2\sqrt{d_F})\epsilon$. Let $x_1, \dots, x_{d_F}$ be an orthonormal basis of $F$ so that, for any $w \in E$, we have

$$ \begin{aligned} \|\ell^T(w)\|_F &= \left( \sum_{k=1}^{d_F} \left| \langle \ell^T(w), x_k \rangle_F \right|^2 \right)^{1/2} \\ &\le \sqrt{d_F} \max_{k,s} \left( s \langle \ell^T(w), x_k \rangle_F \right) \\ &= \max_{k,s} \left\langle w, \ell(s \sqrt{d_F} x_k) \right\rangle_E, \quad (118) \end{aligned} $$

where, here and in the following, $\max_{k,s}$ is the maximum over $s \in \{\pm 1\}$ and $k = 1 \dots d_F$. By definition of $k^{(o)}$ in Definition 1, we also have: for all $\xi' \in D_x(\xi, \epsilon_0)$,

$$ \begin{aligned} e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} k^{(o)}(\ell(\xi'), y) &= h(y) e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} + \langle S(y), \ell(\xi') \rangle_E - A \circ \ell(\xi') \\ &= h(y) e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} + \langle S(y), \ell(\xi' - \xi) \rangle_E + \langle S(y), \ell(\xi) \rangle_E - A \circ \ell(\xi') \end{aligned} $$

Now using that

$$ |\langle S(y), \ell(\xi' - \xi) \rangle_E| = |\langle \ell^T \circ S(y), \xi' - \xi \rangle_F| \le \epsilon_0 \|\ell^T \circ S(y)\|_F $$
(since $\|x\|_F = 1$ and $\xi' \in D_x(\xi, \epsilon_0)$) paired up with (118), we can deduce that

$$
\begin{align}
e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} k^{(o)}(\ell(\xi'), y) &\le h(y) e^{2\epsilon_0 \|\ell^T \circ S(y)\|_F + \langle S(y), \ell(\xi) \rangle_E - A \circ \ell(\xi')} \nonumber \\
&\le \max_{k,s} h(y) \exp \left( \left\langle S(y), \ell(\xi + 2s\sqrt{d_F}\epsilon_0 x_k) \right\rangle_E - A \circ \ell(\xi') \right) \nonumber \\
&\le \max_{k,s} k^{(o)}(\ell(\xi''_{s,k}), y) e^{|A \circ \ell(\xi'') - A \circ \ell(\xi')|}, \tag{119}
\end{align}
$$

where $\xi''_{s,k} = \xi + 2s\sqrt{d_F}\epsilon_0 x_k$. Note that we took $\epsilon$ small enough so that for all $s, k$,

$$
D_x(\xi, \epsilon' + \epsilon) \cup D_{s x_k}(\xi, \epsilon' + \epsilon) \subset B_F \left( \xi, 2 \left( 1 + \sqrt{d_F} \right) \epsilon \right) \subset \Upsilon(\tilde{T}) .
$$

$\xi'$ and $\xi''_{s,k}$ are in $D_x(\xi, \epsilon_0)$ and $D_{s x_k}(\xi, 2\sqrt{d_F}\epsilon_0)$, respectively. Using (116) of Assertion (i) twice with $\epsilon \le \epsilon' < \epsilon''$ successively replaced by $2\sqrt{d_F}\epsilon_0 \le \epsilon' \le 2(1+\sqrt{d_F})\epsilon$ and $\epsilon_0 \le \epsilon' \le 2(1+\sqrt{d_F})\epsilon$, we get that

$$
\begin{align*}
|A \circ \ell (\xi''_{s,k}) - A \circ \ell (\xi')| &\le |A \circ \ell (\xi''_{s,k}) - A \circ \ell (\xi)| + |A \circ \ell (\xi) - A \circ \ell (\xi')| \\
&\le 2\sqrt{d_F} \epsilon_0 \max (\|\ell^\nabla(\xi)\|_F, \|\ell^\nabla(\xi + \epsilon' s x_k)\|_F) \\
&\quad + \epsilon_0 \max (\|\ell^\nabla(\xi)\|_F, \|\ell^\nabla(\xi + \epsilon' x)\|_F) \\
&\le \epsilon_0' C_{\xi,x}(\ell),
\end{align*}
$$

where the last line follows by setting

$$
C_{\xi,x}(\ell) = \max \left\{ \| \ell^{\nabla}(\xi + \epsilon' x') \|_F : x' \in \{0, x, s x_k\}, \text{ with } s \in \{\pm 1\}, k = 1 \dots d_F \right\}.
$$

Combining this inequality with (119) then yields : for all $\xi' \in D_x(\xi, \epsilon_0)$,

$$
\begin{align*}
e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} k^{(o)}(\ell(\xi'), y) &\le e^{\epsilon_0' C_{\xi,x}(\ell)} \max_{k,s} k^{(o)}(\ell(\xi''_{s,k}), y) \\
&\le e^{\epsilon_0' C_{\xi,x}(\ell)} \sum_{k,s} k^{(o)}(\ell(\xi''_{s,k}), y).
\end{align*}
$$

Observe that the right-hand side of the previous display does not depend on $\xi'$ and that the
summands in the sum (which has $2d_F$ terms) all have integrals equal w.r.t. $\nu$ equal to 1.
The bound (117) follows and Assertion (ii) is proved.
$\blacksquare$

**Lemma 12** Consider an ELM family as in Definition 2. Assume (C1) and (C3). Let $\vartheta_0 \in \tilde{T}$, $\tau_0 \in M$ and $\tilde{\varphi}: \mathcal{Y} \to \mathbb{R}_+$ be a density function such that

$$
\int_{\mathcal{Y}} \tilde{\varphi}(y) e^{\tilde{s}_{\vartheta_0, \tau_0}(y)} \nu(dy) < \infty. \quad (120)
$$

Define $\varphi(\ell, y)$, $\tilde{\varphi}(\ell)$ and $w_{\tau_0}^\nabla$ as in (52), (59) and (60). Then for all $\xi \in \Upsilon(\tilde{T})$ and $x \in F$ such that $\|x\|_F = 1$, there exists $\epsilon_0 > 0$ such that

$$
\int_{\mathcal{Y} \times \mathbb{Z}} \| \ell^T \circ S(y) \|_F \sup_{\xi' \in D_x(\xi, \epsilon_0)} \left( k^{(o)}(\ell(\xi'), y) \right) \tilde{\varphi}(\ell) \nu(dy) \tau_0(d\ell) < \infty. \quad (121)
$$
Moreover, $w_{\tau_0}^{\nabla}$ is continuous on $\Upsilon(\tilde{T})$ in every direction, and, for all $\vartheta \in \tilde{T}$, it holds that

$$ w_{\tau_0}^{\nabla} \circ \Upsilon(\vartheta) = \mathbb{E}_{\vartheta, \tau_0} [L^T \circ S(Y) \tilde{\varphi}(L)] . \qquad (122) $$

**Proof** Let $\ell \in \mathbb{Z}$ is such that $\ell \circ \Upsilon(\tilde{T}) \subset \text{Int} E_0$, $\xi \in \Upsilon(\tilde{T})$ and $x \in F$ such that $\|x\|_F = 1$. Then for $\epsilon > 0$ small enough, we have $B_x(\xi, 2(1+\sqrt{d_F})\epsilon) \subset \Upsilon(\tilde{T})$, and for all $\epsilon_0 \in (0, \epsilon]$,

$$ \int_{\Upsilon} e^{\epsilon_0 \|\ell^T \circ S(y)\|_F} \sup_{\xi' \in D_x(\xi, \epsilon_0)} \left( k^{(o)}(\ell(\xi'), y) \right) \nu(dy) \le 2 d_F \sum_{\xi' \in \Delta} e^{\epsilon_0' \|\ell^{\nabla}(\xi')\|_F}, $$

where we used Assertion (ii) of Lemma 11 with $\epsilon' = (1 + 2\sqrt{d_F})\epsilon$, $\epsilon_0' = (1 + 2\sqrt{d_F})\epsilon_0$ and $\Delta$ is a set of $2(d_F + 1)$ points in $B_x(\xi, 2(1+\sqrt{d_F})\epsilon)$ only depending on $\xi$, $x$ and $\epsilon$. Since $\|\ell^T \circ S(y)\|_F \le \frac{1}{\epsilon_0} e^{\epsilon_0 \|\ell^T \circ S(y)\|_F}$, we get that:

$$ \begin{aligned} G_{\xi, \epsilon_0}(\ell) &:= \int_{\Upsilon} \|\ell^T \circ S(y)\|_F \sup_{\xi' \in D_x(\xi, \epsilon_0)} \left( k^{(o)}(\ell(\xi'), y) \right) \nu(dy) \\ &\le \frac{2 d_F}{\epsilon_0} \sum_{\xi' \in \Delta} e^{\epsilon_0' \|\ell^{\nabla}(\xi')\|_F}. \end{aligned} $$

In particular by (C1), with Lemma 9, this bounds hold for $\tau_0$-almost all $\ell \in \mathbb{Z}$. By (C3), since $\Delta$ is a finite set only depending on $\xi$, $x$ and $\epsilon$ and satisfying $\Delta \subset B_x(\xi, 2(1+\sqrt{d_F})\epsilon) \subset \Upsilon(\tilde{T})$, there exist $\epsilon_\Delta$, $C_\Delta > 0$ only depending on $\xi$, $x$ and $\epsilon$ such that, for all $\xi' \in \Delta$,

$$ \mathbb{E}_{\vartheta_0, \tau_0} \left[ e^{\epsilon_\Delta \|L^{\nabla}(\xi')\|_F} \mid Y \right] \le C_\Delta e^{\tilde{s}_{\vartheta_0, \tau_0}(Y)} \mathbb{P}_{\vartheta_0, \tau_0} - \text{a.s.} $$

Since $\Delta$ has cardinal $2(d_F + 1)$, taking $\epsilon_0 \in (0, \epsilon]$ small enough to have $\epsilon_0' \le \epsilon_\Delta$, this gives us that

$$ g_{\xi, \epsilon_0}(Y) := \mathbb{E}_{\vartheta_0, \tau_0} [G_{\xi, \epsilon_0}(L) | Y] \le \frac{4 C_\Delta d_F(d_F + 1)}{\epsilon_0} e^{\tilde{s}_{\vartheta_0, \tau_0}(Y)} \mathbb{P}_{\vartheta_0, \tau_0} - \text{a.s.} $$

Now observing that the ratio in (52) is the conditional density of $L$ applied to $\ell$ given $Y = y$ under $\mathbb{P}_{\vartheta_0, \tau_0}$ and that the marginal distribution of $\mathbb{P}_{\vartheta_0, \tau_0}$ on $\Upsilon$ is equivalent to $\nu$, we obtain that

$$ \begin{aligned} \int_{\mathbb{Z}} G_{\xi, \epsilon_0}(\ell) \tilde{\varphi}(\ell) \tau_0(d\ell) &= \int_{\mathbb{Z} \times \Upsilon} G_{\xi, \epsilon_0}(\ell) \varphi(\ell, y) \tau_0(d\ell) \nu(dy) \\ &= \int_{\Upsilon} g_{\xi, \epsilon_0}(y) \tilde{\varphi}(y) \nu(dy) \\ &\le \frac{4 C_\Delta d_F(d_F + 1)}{\epsilon_0} \int_{\Upsilon} e^{\tilde{s}_{\vartheta_0, \tau_0}(y)} \tilde{\varphi}(y) \nu(dy), \end{aligned} $$

Hence, using condition (120), we deduce that

$$ \int_{\mathbb{Z}} G_{\xi, \epsilon_0}(\ell) \tilde{\varphi}(\ell) \tau_0(d\ell) < \infty, $$