# 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, United Kingdom
<EMAIL_ADDRESS>
&Randal Douc
SAMOVAR, Télécom SudParis
Institut Polytechnique de Paris, France
<EMAIL_ADDRESS>
François Roueff
LTCI, Télécom Paris
Institut Polytechnique de Paris, France
<EMAIL_ADDRESS>

###### Abstract

In this paper, we introduce a novel family of iterative algorithms which carry
out $\alpha$-divergence minimisation in a Variational Inference context. They
do so by ensuring a systematic decrease at each step in the
$\alpha$-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. Notably, our approach permits to build on
various methods previously proposed for $\alpha$-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.

## 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 1, 2] have attracted a lot of attention as they have empirically been
shown to be widely applicable to many high-dimensional machine-learning
problems [3, 4, 5].

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 [6, 7].

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 [8, 9, 10]. This is especially inconvenient when the
variational family $\mathcal{Q}$ does not exactly match the posterior density
[11, 12] and even more so if one wants to appeal to Importance Sampling
methods to approximate integrals of interest in a Bayesian Inference setting
[9, 10].

To avoid this hurdle, advances in Variational Inference turned to more general
families of divergences such as the $\alpha$-divergence [13, 14] and Rényi’s
$\alpha$-divergence [15, 16]. 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$ [8]. For this reason, they have notably been used in [17,
18, 19, 20, 21, 22].

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:

  * $\bullet$

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

  * $\bullet$

In Section 3, we consider the typical Variational Inference case where $q$
belongs to a parametric family. In this particular case, we state in 1
conditions which ensure a systematic decrease in the $\alpha$-divergence at
each step for all $\alpha\in[0,1)$. We then show in 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 1 and 1. Furthermore, we derive in 2 additional
iterative schemes satisfying the conditions of 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.

  * $\bullet$

In Section 4, we further 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 both 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 4 and 5, which respectively generalise the maximisation approach of 1
and the gradient-based approach of 2 to the case of mixture models.

  * $\bullet$

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. Notably, we connect our approach to the
Power Descent algorithm from [22] and provide in 2 additional monotonicity
results which go beyond the case $\alpha\in[0,1)$. We also obtain that the
Mixture Population Monte Carlo (M-PMC) algorithm, an integrated Expectation
Maximization algorithm introduced in [23], can be recovered as a special case
of our framework.

  * $\bullet$

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 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 5).

  * $\bullet$

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

  * $\bullet$

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.

## 2 Notation and optimisation problem

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

$p(y|\mathscr{D})=\frac{p(y,\mathscr{D})}{p(\mathscr{D})}=\frac{p_{0}(y)p(\mathscr{D}|y)}{p(\mathscr{D})}\;,$

where the normalisation constant
$p(\mathscr{D})=\int_{\mathsf{Y}}p_{0}(y)p(\mathscr{D}|y)\nu(\mathrm{d}y)$ 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 8].

More precisely, let us denote by $\mathbb{P}$ the probability measure on
$(\mathsf{Y},\mathcal{Y})$ with corresponding density $p(\cdot|\mathscr{D})$
with respect to $\nu$. As for the variational family, let us denote by
$\mathbb{Q}$ the probability measure on $(\mathsf{Y},\mathcal{Y})$ with
associated density $q\in\mathcal{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)}\left[u^{\alpha}-1\right]$ for all
$\alpha\in\mathbb{R}\setminus\left\\{0,1\right\\}$. 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 24]
is given by

$\displaystyle
D_{\alpha}(\mathbb{Q}||\mathbb{P})=\int_{\mathsf{Y}}f_{\alpha}\left(\frac{q(y)}{p(y|\mathscr{D})}\right)p(y|\mathscr{D})\nu(\mathrm{d}y)\;,$
(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 the optimisation problem above is
equivalent to solving

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

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

$\displaystyle\Psi_{\alpha}(q;p)=\int_{\mathsf{Y}}f_{\alpha}\left(\frac{q(y)}{p(y)}\right)p(y)\nu(\mathrm{d}y)\;.$

As the normalisation constant does not appear anymore in the optimisation
problem (2), this formulation is often preferred in practice. Therefore, we
consider the general optimisation problem

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

where $p$ is any measurable positive function on $(\mathsf{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 $({\mathsf{T}},\mathcal{T})$ be a
measurable space, $K:(\theta,A)\mapsto\int_{A}k(\theta,y)\nu(\mathrm{d}y)$ be
a Markov transition kernel on ${\mathsf{T}}\times\mathcal{Y}$ with kernel
density $k$ defined on ${\mathsf{T}}\times\mathsf{Y}$, we consider a
parametric family of the form

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

## 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.

  1. (A1)

The kernel density $k$ on ${\mathsf{T}}\times\mathsf{Y}$, the function $p$ on
$\mathsf{Y}$ and the $\sigma$-finite measure $\nu$ on
$(\mathsf{Y},\mathcal{Y})$ satisfy, for all
$(\theta,y)\in{\mathsf{T}}\times\mathsf{Y}$, $k(\theta,y)>0$, $p(y)\geqslant
0$ and $0<\int_{\mathsf{Y}}p(y)\nu(\mathrm{d}y)<\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\geqslant 1}$ valued in
${\mathsf{T}}$ such that the sequence
$(\Psi_{\alpha}(k(\theta_{n},\cdot))_{n\geqslant 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^{\prime}\in{\mathsf{T}}$ such that
$\Psi_{\alpha}(k(\theta^{\prime},\cdot))<\infty$, it holds that

$\displaystyle\Psi_{\alpha}(k(\theta,\cdot))$
$\displaystyle\leqslant\int_{\mathsf{Y}}\frac{k(\theta^{\prime},y)^{\alpha}p(y)^{1-\alpha}}{\alpha-1}\log\left(\frac{k(\theta,y)}{k(\theta^{\prime},y)}\right)\nu(\mathrm{d}y)+\Psi_{\alpha}(k(\theta^{\prime},\cdot))\;.$
(4)

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

###### Proof.

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

  1. (a)

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

$\displaystyle\Psi_{0}(k(\theta,\cdot))=-\int_{\mathsf{Y}}p(y)\log\left(\frac{k(\theta,y)}{k(\theta^{\prime},y)}\right)\nu(\mathrm{d}y)+\Psi_{0}(k(\theta^{\prime},\cdot))\;.$

  2. (b)

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

$\displaystyle\Psi_{\alpha}(k(\theta,\cdot))$
$\displaystyle=\int_{\mathsf{Y}}\frac{\left[\left(\frac{k(\theta,y)}{p(y)}\right)^{\alpha}-1\right]}{\alpha(\alpha-1)}p(y)\nu(\mathrm{d}y)$
$\displaystyle=\int_{\mathsf{Y}}\left(\frac{k(\theta^{\prime},y)}{p(y)}\right)^{\alpha}\frac{\left[\left(\frac{k(\theta,y)}{k(\theta^{\prime},y)}\right)^{\alpha}-1\right]}{\alpha(\alpha-1)}p(y)\nu(\mathrm{d}y)+\Psi_{\alpha}(k(\theta^{\prime},\cdot))$

Furthermore, using that $\log(u^{\alpha})\leqslant 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)}\leqslant\frac{\log(u^{\alpha})}{\alpha(\alpha-1)}=\frac{\log(u)}{\alpha-1}\;.$

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

$\displaystyle\Psi_{\alpha}(k(\theta,\cdot))$
$\displaystyle\leqslant\int_{\mathsf{Y}}\frac{k(\theta^{\prime},y)^{\alpha}p(y)^{1-\alpha}}{\alpha-1}\log\left(\frac{k(\theta,y)}{k(\theta^{\prime},y)}\right)\nu(\mathrm{d}y)+\Psi_{\alpha}(k(\theta^{\prime},\cdot))$

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

∎

This result then allows us to deduce 1 below.

######  Theorem 1.

Assume (A1). Let $\alpha\in[0,1)$ and starting from an initial
$\theta_{1}\in{\mathsf{T}}$, let $(\theta_{n})_{n\geqslant 1}$ be defined
iteratively such that for all $n\geqslant 1$,

$\displaystyle\int_{\mathsf{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(\mathrm{d}y)\leqslant
0\;.$ (5)

Further assume that $\Psi_{\alpha}(k(\theta_{1},\cdot))<\infty$. Then,
$\Psi_{\alpha}(k(\theta_{n+1},\cdot))\leqslant\Psi_{\alpha}(k(\theta_{n},\cdot))$
for all $n\geqslant 1$ with an equality occurring only if (5) is also an
equality.

###### Proof.

Applying 1 iteratively for $n=1,2,\dots$ with $\theta=\theta_{n+1}$ and
$\theta^{\prime}=\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))\leqslant\int_{\mathsf{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(\mathrm{d}y)+\Psi_{\alpha}(k(\theta_{n},\cdot))\leqslant\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\geqslant 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^{\prime}=\theta_{n}$. ∎

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

$\displaystyle\ratio(y)={k(\theta_{n},y)^{\alpha}p(y)^{1-\alpha}}$ (6)
$\displaystyle\check{\varphi}_{n}^{(\alpha)}=\frac{\ratio}{\int\ratio\;\mathrm{d}\nu}\;.$

Note that under (A1), the normalisation in $\check{\varphi}_{n}^{(\alpha)}$
satisfies $0<\int\ratio\;\mathrm{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\geqslant 1}$ be a non-
negative sequence. Starting from an initial $\theta_{1}\in{\mathsf{T}}$, let
$(\theta_{n})_{n\geqslant 1}$ be defined iteratively as follows

$\displaystyle\theta_{n+1}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{Y}}\left[\ratio(y)+b_{n}k(\theta_{n},y)\right]\log\left(\frac{k(\theta,y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\;,\quad
n\geqslant 1\;{{,}}$ (7)

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

###### Proof.

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

$\displaystyle\int_{\mathsf{Y}}\ratio(y)\log\left(\frac{k(\theta_{n+1},y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\geqslant
0\;.$

Since
$\int_{\mathsf{Y}}k(\theta_{n},y)\log\left(\frac{k(\theta_{n+1},y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)=-D_{1}(K(\theta_{n},\cdot)||K(\theta_{n+1},\cdot))\leqslant
0$, we also obtain that

$\displaystyle\int_{\mathsf{Y}}\ratio(y)\log\left(\frac{k(\theta_{n+1},y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\geqslant\int_{\mathsf{Y}}\left[\ratio(y)+b_{n}k(\theta_{n},y)\right]\log\left(\frac{k(\theta_{n+1},y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\;.$

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

$\displaystyle\int_{\mathsf{Y}}\left[\ratio(y)+b_{n}k(\theta_{n},y)\right]\log\left(\frac{k(\theta_{n+1},y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\geqslant
0\;,\quad n\geqslant 1\;,$

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

We now make three comments regarding 1.

  * •

Firstly, the r.h.s of (7) can be simplified, yielding the equivalent iterative
scheme

$\displaystyle\theta_{n+1}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{Y}}\left[\ratio(y)+b_{n}k(\theta_{n},y)\right]\log
k(\theta,y)\nu(\mathrm{d}y)\;,\quad n\geqslant 1\;,$

provided that
$\int_{\mathsf{Y}}\left[\ratio(y)+b_{n}k(\theta_{n},y)\right]|\log
k(\theta_{n},y)|\nu(\mathrm{d}y)<\infty$ for all $n\geqslant 1$.

  * •

Secondly, a remarkable aspect is 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\\}$. This will allow us to solve the argmax problem in
(7) for some notable choices of kernel density $k$ and thus to use the
maximisation approach of 1 to derive simple update rules for
$(\theta_{n})_{n\geqslant 1}$. An example and a lemma are provided thereafter
to illustrate this fact and more general results will follow later in Section
6.

  * •

Thirdly, the sequence $(b_{n})_{n\geqslant 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 provide below a first motivating example where the maximisation approach of
1 is applicable.

######  Example 1 (Maximisation approach for a Gaussian density).

We consider a $d$-dimen-sional Gaussian density, in which case
$\mathsf{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{\mathsf{T}}$ denotes the mean and covariance matrix of
the Gaussian density. We let the parameter set ${\mathsf{T}}$ include all
possible means $m$ in $\mathbb{R}^{d}$ and all possible positive-definite
covariance matrices $\Sigma$. Finally, denoting by $\left\|\cdot\right\|$ the
Euclidian norm, we assume that the non-negative function $p$ defined on
$\mathsf{Y}$ satisfies

$\displaystyle
0<\int_{\mathsf{Y}}\left(1+\|y\|^{2/(1-\alpha)}\right)\,p(y)\;\mathrm{d}y<\infty\;.$
(8)

Then (A1) holds. Furthermore, starting from any
$\theta_{1}=(m_{1},\Sigma_{1})\in{\mathsf{T}}$ so that
$\Psi_{\alpha}(k(\theta_{1},\cdot))<\infty$ and using the notation
$\theta_{n}=(m_{n},\Sigma_{n})$ for all $n\geqslant 1$, it holds that: for all
$n\geqslant 1$ and all $\gamma_{n}\in(0,1]$, there exists $b_{n}\geqslant 0$
such that the argmax problem (7) has a unique solution defined by

$\displaystyle\begin{split}&m_{n+1}=\gamma_{n}\,\hat{m}_{n}+(1-\gamma_{n})\,m_{n}\;,\\\
&\Sigma_{n+1}=\gamma_{n}\,\widehat{\Sigma}_{n}+(1-\gamma_{n})\,\Sigma_{n}+\gamma_{n}(1-\gamma_{n})\,\left(\widehat{m}_{n}-m_{n}\right)\left(\widehat{m}_{n}-m_{n}\right)^{T}\;,\end{split}$
(9)

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

$\displaystyle\widehat{m}_{n}=\int_{\mathsf{Y}}y\,\check{\varphi}_{n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)\;,$
$\displaystyle\widehat{\Sigma}_{n}=\int_{\mathsf{Y}}yy^{T}\,\check{\varphi}_{n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)-\widehat{m}_{n}\widehat{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 1 (in particular we will see in 6 that
the case $\gamma_{n}=1$ corresponding to $b_{n}=0$ will require an additional
non-degenerate condition on $p$.)

Observe then that, by 1, the systematic decrease in
$(\Psi_{\alpha}(k(\theta_{n},\cdot)))_{n\geqslant 1}$ holds for any choice of
sequence $(\gamma_{n})_{n\geqslant 1}$ valued in $(0,1]$ in the update
formulas (9). In addition, $\gamma_{n}$ permits us to build a tradeoff at time
$n\geqslant 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 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. For this choice of
variational family, the following lemma indeed states that the global argmax
problem can be separated into component-wise argmax problems.

######  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)},\ldots,\theta^{(L)})\in{\mathsf{T}}$ and
$y=(y^{(1)},\ldots,y^{(L)})\in\mathsf{Y}$. Then, starting from
$\theta_{1}\in{\mathsf{T}}$ and using the notation
$\theta_{n}=(\theta_{n}^{(1)},\ldots,\theta_{n}^{(L)})$ for all $n\geqslant
1$, solving (7) yields the following update formulas: for all $n\geqslant 1$,

$\displaystyle\theta^{(\ell)}_{n+1}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta^{(\ell)}}\int_{\mathsf{Y}}\left[\ratio(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(\mathrm{d}y)\;,\quad{{\ell=1\ldots
L}}\;.$

Interestingly, it is also possible to derive iterative schemes satisfying (5)
under additional smoothness conditions (see Section A.1 for the definition of
$\beta$-smoothness), as written in 2.

######  Corollary 2 (Gradient-based approach).

Assume (A1). Let ${\mathsf{T}}\subseteq\mathbb{R}^{d}$ be a convex set, let
$\alpha\in[0,1)$ and let $(\gamma_{n})_{n\geqslant 1}$ be valued in $(0,1]$.
Starting from an initial $\theta_{1}\in{\mathsf{T}}$, let
$(\theta_{n})_{n\geqslant 1}$ be defined iteratively as follows

$\displaystyle\theta_{n+1}$
$\displaystyle=\theta_{n}-\frac{\gamma_{n}}{\beta_{n}}\nabla
g_{n}(\theta_{n})\;,\quad n\geqslant 1\;,$ (10)

where $(g_{n})_{n\geqslant 1}$ is the sequence of functions defined by: for
all $n\geqslant 1$ and all $\theta\in{\mathsf{T}}$

$\displaystyle
g_{n}(\theta)=\int_{\mathsf{Y}}\frac{{{\ratio(y)}}}{\alpha-1}\log\left(\frac{k(\theta,y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\;,$
(11)

and $g_{n}$ is assumed to be $\beta_{n}$-smooth. Then (5) holds for all
$n\geqslant 1$ and we can apply 1.

###### Proof.

Since $\gamma_{n}\in(0,1]$ and $g_{n}$ is a $\beta_{n}$-smooth function by
assumption, we can apply 2 and we obtain that for all $n\geqslant 1$,

$\displaystyle
g_{n}(\theta_{n})-g_{n}\left(\theta_{n}-\frac{\gamma_{n}}{\beta_{n}}\nabla
g_{n}(\theta_{n})\right)\geqslant\frac{\gamma_{n}}{2\beta_{n}}\left\|\nabla
g_{n}(\theta_{n})\right\|^{2}\;.$

Thus, by definition of $\theta_{n+1}$ in (10), we have

$0=g_{n}(\theta_{n})\geqslant g_{n}(\theta_{n+1})\;,$

which in turn implies (5) and the proof is concluded. ∎

Let us now reflect on the implications of 2. Under common differentiability
assumptions, we can write: for all $n\geqslant 1$ and all
${{\theta^{\prime}}}\in{\mathsf{T}}$

$\displaystyle\nabla
g_{n}(\theta^{\prime})=\int_{\mathsf{Y}}\frac{{{\ratio(y)}}}{\alpha-1}\left.\frac{\partial\log
k(\theta,y)}{\partial\theta}\right|_{(\theta,y)=(\theta^{\prime},y)}\nu(\mathrm{d}y)\;.$

Then, (10) becomes

$\displaystyle\theta_{n+1}$
$\displaystyle=\theta_{n}-\frac{\gamma_{n}}{\beta_{n}}\int_{\mathsf{Y}}\frac{{{\ratio(y)}}}{\alpha-1}\left.\frac{\partial\log
k(\theta,y)}{\partial\theta}\right|_{(\theta,y)=(\theta_{n},y)}\nu(\mathrm{d}y)\;,\quad
n\geqslant 1\;.$ (12)

Written under this form, one may observe that the iterative scheme (12) bears
similarities with Gradient Descent iterations for $\alpha$-divergence and
Rényi’s $\alpha$-divergence minimisation.

Indeed, given a well-chosen learning rate policy $(r_{n})_{n\geqslant 1}$ and
setting $p=p(\cdot,\mathscr{D})$, such Gradient Descent iterations are given
respectively by

$\displaystyle\theta_{n+1}=\theta_{n}-r_{n}\int_{\mathsf{Y}}\frac{\ratio(y)}{\alpha-1}\left.\frac{\partial\log
k(\theta,y)}{\partial\theta}\right|_{(\theta,y)=(\theta_{n},y)}\nu(\mathrm{d}y)\;,\quad
n\geqslant 1$
$\displaystyle\theta_{n+1}=\theta_{n}-r_{n}\int_{\mathsf{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(\mathrm{d}y)\;,\quad
n\geqslant 1$

(we refer to Section 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\geqslant 1}$ from 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)=\mathcal{N}(y;m,\Sigma)$, where this time $\theta=m$. Further
assume that ${\mathsf{T}}\subseteq\mathbb{R}^{d}$ is a convex subset and that
$\Sigma=\sigma^{2}\boldsymbol{I_{d}}$, where $\sigma^{2}>0$ and
$\boldsymbol{I_{d}}$ denotes the $d$-dimensional identity matrix. Finally,
denoting by $\left\|\cdot\right\|$ the Euclidian norm, we assume that the non-
negative function $p$ defined on $\mathsf{Y}$ satisfies

$\displaystyle
0<\int_{\mathsf{Y}}\left(1+\|y\|^{1/(1-\alpha)}\right)p(y)\mathrm{d}y<\infty\;.$
(13)

Then, (A1) holds. Furthermore, setting
$\beta_{n}=\sigma^{-2}(1-\alpha)^{-1}\int\ratio\mathrm{d}\nu$, the conditions
on $(g_{n})_{n\geqslant 1}$ from 2 are satisfied and using the notation
$\theta_{n}=m_{n}$, the mean of the Gaussian density $\mathcal{N}(y;m,\Sigma)$
can be optimised using the iterative scheme given by: for all $n\geqslant 1$
and all $\gamma_{n}\in[0,1)$

$\displaystyle{{m_{n+1}}}$
$\displaystyle=\gamma_{n}\hat{m}_{n}+\left(1-\gamma_{n}\right)m_{n}\;,$ (14)

where $\hat{m}_{n}=\int_{\mathsf{Y}}y\leavevmode\nobreak\
\check{\varphi}_{n}^{(\alpha)}(y)\leavevmode\nobreak\ \nu(\mathrm{d}y)$ (we
refer to 8 for the detailed derivation). Setting $p=p(\cdot,\mathscr{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.

  * •

2 sheds light on the links between our approach and the more traditional
gradient descent methodology for optimising objective functions in Variational
Inference [19, 20].

Unlike the usual gradient descent methodology, 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 25, 26]. Yet, it can be difficult to satisfy those smoothness
asumptions in practice, even if some results have been derived for specific
variational families when using the exclusive Kullback-Leibler divergence as
the objective function [27].

To the best of our knowledge, no such results have been proven for
$\alpha$-divergence minimisation. In our case, we obtain that a smoothness
condition on $g_{n}$ translates into a systematic decrease in $\Psi_{\alpha}$.
As we shall see in the forthcoming section, being able to establish
monotonicity results on $(\theta_{n})_{n\geqslant 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 1 considerably eased the derivation
of iterative update formulas for well-chosen kernel densities $k$ that can
exploit this log structure (1-2 and 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 fomulas 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 1 to Gaussian Mixture Models).

As we shall see in the next section, further theoretical developpements 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}^{\star}$, we introduce the simplex of
$\mathbb{R}^{J}$:

$\mathcal{S}_{J}=\left\\{{\boldsymbol{\lambda}}=(\lambda_{1},\ldots,\lambda_{J})\in\mathbb{R}^{J}\;:\;\forall
j\in\left\\{1,\ldots,J\right\\},\;\lambda_{j}\geqslant
0\;\mbox{and}\;\sum_{j=1}^{J}\lambda_{j}=1\right\\}\;,$

we define
$\mathcal{S}_{J}^{+}=\left\\{{\boldsymbol{\lambda}}\in\mathcal{S}_{J}\;:\;\forall
j\in\left\\{1,\ldots,J\right\\},\;\lambda_{j}>0\right\\}$ and we denote
$\Theta=(\theta_{1},\ldots,\theta_{J})\in{\mathsf{T}}^{J}$. We consider the
mixture model variational family given by

$\displaystyle\mathcal{Q}=\left\\{q:y\mapsto\mu_{{\boldsymbol{\lambda}},\Theta}k(y):=\sum_{j=1}^{J}\lambda_{j}k(\theta_{j},y)\;:\;{\boldsymbol{\lambda}}\in\mathcal{S}_{J},\Theta\in{\mathsf{T}}^{J}\right\\}$
(15)

that is, we are interested in solving the optimisation problem

$\displaystyle\inf_{{\boldsymbol{\lambda}}\in\mathcal{S}_{J},\Theta\in{\mathsf{T}}^{J}}\Psi_{\alpha}(\mu_{{\boldsymbol{\lambda}},\Theta}k)\;,$

with $J>1$ (notice in particular that the framework from Section 3 would have
corresponded to having taken $J=1$ in the optimisation problem above). Let us
next denote ${\boldsymbol{\lambda}_{n}}=(\lambda_{j,n})_{1\leqslant j\leqslant
J}$ and $\Theta_{n}=(\theta_{j,n})_{1\leqslant j\leqslant J}$ for all
$n\geqslant 1$. For convenience, we also introduce the shorthand notation
$\mu_{n}k=\mu_{{\boldsymbol{\lambda}_{n}},\Theta_{n}}k$ and

$\displaystyle\varphi_{j,n}^{(\alpha)}(y)=k(\theta_{j,n},y)\left(\frac{\mu_{n}k(y)}{p(y)}\right)^{\alpha-1}$
(16)
$\displaystyle\check{\varphi}_{j,n}^{(\alpha)}=\frac{\varphi_{j,n}^{(\alpha)}}{\int\varphi_{j,n}^{(\alpha)}\mathrm{d}\nu}$

for all $\alpha\in[0,1)$, all $j=1\ldots J$, all $n\geqslant 1$ and all
$y\in\mathsf{Y}$.

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

######  Corollary 3.

Assume (A1). Let $J\in\mathbb{N}^{\star}$, let $\alpha\in[0,1)$ and starting
from an initial parameter set
$({\boldsymbol{\lambda}_{1}},\Theta_{1})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$,
let $({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ be defined
iteratively such that for all $n\geqslant 1$,

$\displaystyle\int_{\mathsf{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(\mathrm{d}y)\leqslant
0\;.$ (17)

Further assume that $\Psi_{\alpha}(\mu_{1}k)<\infty$. Then,
$\Psi_{\alpha}(\mu_{n+1}k)\leqslant\Psi_{\alpha}(\mu_{n}k)$ for all
$n\geqslant 1$.

Observe that we are in a less favourable situation in 3 with $J>1$ compared to
the case $J=1$ we studied previously in Section 3. Indeed, we now have a ratio
of sums inside the $\log$ function in (17), meaning that the approach from 1
to derive simple iterative schemes does not directly transfer to the
variational family (15) for choices of kernel densities $k$ identified in 1-2
and 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}^{\star}$, $\alpha\in[0,1)$ and starting from
an initial parameter set
$({\boldsymbol{\lambda}_{1}},\Theta_{1})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$,
let $({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ be defined
iteratively such that for all $n\geqslant 1$,

$\displaystyle\int_{\mathsf{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(\mathrm{d}y)\leqslant
0$ (18)
$\displaystyle\int_{\mathsf{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(\mathrm{d}y)\leqslant
0\;.$ (19)

Further assume that $\Psi_{\alpha}(\mu_{1}k)<\infty$. Then,
$\Psi_{\alpha}(\mu_{n+1}k)\leqslant\Psi_{\alpha}(\mu_{n}k)$ for all
$n\geqslant 1$.

###### Proof.

By 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\mathsf{Y}$ and all $n\geqslant 1$,

$\displaystyle\frac{1}{\alpha-1}\log\left(\frac{\mu_{n+1}k(y)}{\mu_{n}k(y)}\right)$
$\displaystyle=\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)$
$\displaystyle\leqslant\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)$
(20)

that is:

$\displaystyle\frac{1}{\alpha-1}\log\left(\frac{\mu_{n+1}k(y)}{\mu_{n}k(y)}\right)\leqslant\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(\mathrm{d}y)$ and using the definition of
$\varphi_{j,n}^{(\alpha)}(y)$ in (16), this in turn implies that: for all
$n\geqslant 1$,

$\int_{\mathsf{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(\mathrm{d}y)\\\
\leqslant\int_{\mathsf{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(\mathrm{d}y)\;.$
(21)

As a consequence, the condition

$\displaystyle\int_{\mathsf{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(\mathrm{d}y)\leqslant
0$ (22)

implies (17). The condition (22) 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
${\boldsymbol{\lambda}_{n+1}}$ in 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
laterly in Section 8, will lead to reduced computational power in practice.

In addition, we have recovered in (19) the key property from the case $J=1$ in
Section 3: compared to the initial condition (17) which involved the $\log$ of
weighted sums of kernel densities, (19) now considers a sum of logs of kernel
densities. This suggests that we will be able to 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 will induce numerical advantages later
illustrated in Section 8.

### 4.1 Choice of $({\boldsymbol{\lambda}_{n}})_{n\geqslant 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\geqslant 1}$ be valued
in $(0,1]$ and let $(\kappa_{n})_{n\geqslant 1}$ be such that
$(\alpha-1)\kappa_{n}\geqslant 0$ at all times $n$. Furthermore, let
$(\Theta_{n})_{n\geqslant 1}$ be any sequence valued in ${\mathsf{T}}^{J}$.
Starting from an initial ${\boldsymbol{\lambda}_{1}}\in\mathcal{S}_{J}^{+}$,
let $({\boldsymbol{\lambda}_{n}})_{n\geqslant 1}$ be defined iteratively such
that for all $n\geqslant 1$

$\displaystyle\lambda_{j,n+1}$
$\displaystyle=\frac{\lambda_{j,n}\left[\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}{\sum_{\ell=1}^{J}\lambda_{\ell,n}\left[\int_{\mathsf{Y}}\varphi_{\ell,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}\;,\quad
j=1\ldots J\;.$ (23)

Then (18) holds.

###### Proof.

We first check that the integrals appearing in (23) are finite:

$\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)<\infty\;,\qquad
j=1,\dots,J\,,\;n\geqslant 1\;.$ (24)

Using (16), that $\lambda_{j,n}>0$, that
$\lambda_{j,n}k(\theta_{j,n},y)\leqslant\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\geqslant 1$,

$\displaystyle\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)$
$\displaystyle=\int_{\mathsf{Y}}k(\theta_{j,n},y)\left(\frac{\mu_{n}k(y)}{p(y)}\right)^{\alpha-1}\nu(\mathrm{d}y)$
$\displaystyle\leqslant\frac{1}{\lambda_{j,n}}\int_{\mathsf{Y}}\mu_{n}k(y)\left(\frac{p(y)}{\mu_{n}k(y)}\right)^{1-\alpha}\nu(\mathrm{d}y)$
$\displaystyle\leqslant\frac{1}{\lambda_{j,n}}\left(\int_{\mathsf{Y}}p(y)\nu(\mathrm{d}y)\right)^{1-\alpha}\;.$

The bound (24) follows by (A1). Now, to prove (18), we treat the cases
$\eta_{n}=1$ and $\eta_{n}\in(0,1)$ separately.

  1. (a)

Case $\eta_{n}=1$. Since $(\alpha-1)\kappa_{n}\geqslant 0$ with
$\alpha\in(0,1)$, we have that

$\kappa_{n}\sum_{j=1}^{J}\lambda_{j,n}\log(\lambda_{j}/\lambda_{j,n})\geqslant
0$

where we have used that
$\sum_{j=1}^{J}\lambda_{j,n}\log(\lambda_{j}/\lambda_{j,n})\leqslant\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

$\displaystyle\int_{\mathsf{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(\mathrm{d}y)+\kappa_{n}\sum_{j=1}^{J}\lambda_{j,n}\log\left(\frac{\lambda_{j,n+1}}{\lambda_{j,n}}\right)\leqslant
0$

that is, since (23) holds,

$\displaystyle\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\log\left(\frac{\lambda_{j,n+1}}{\lambda_{j,n}}\right)\leqslant
0\;.$ (25)

Notice then that by definition of $(\lambda_{j,n+1})_{1\leqslant j\leqslant
J}$ when $\eta_{n}=1$, we can write

$\displaystyle{\boldsymbol{\lambda}_{n+1}}=\mathop{\mathrm{argmin}}_{{\boldsymbol{\lambda}}\in\mathcal{S}_{J}^{+}}\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\log\left(\frac{\lambda_{j}}{\lambda_{j,n}}\right)\;.$

[Indeed, setting
$\beta_{j}=\lambda_{j,n}\left[\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]$
and $\bar{\beta}_{j}=\beta_{j}/\sum_{\ell=1}^{J}\beta_{\ell}$ for all
$j=1\ldots J$, we have that
$\sum_{j=1}^{J}\bar{\beta}_{j}\log\left(\bar{\beta}_{j}/\lambda_{j}\right)\geqslant
0$ and that this quantity is minimal when $\lambda_{j}=\bar{\beta}_{j}$ for
$j=1\ldots J$.] This implies (25) and settles the case $\eta_{n}=1$.

  2. (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})\geqslant\frac{1}{\epsilon}\left(1-\frac{1}{u^{\epsilon}}\right)\;.$

Indeed, since
$\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\leqslant
0$ for all $j=1\ldots J$, we can then write that for all $\epsilon>0$,

$\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\log\left(\frac{\lambda_{j,n+1}}{\lambda_{j,n}}\right)\\\
\leqslant\frac{1}{\epsilon}\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\left[1-\left(\frac{\lambda_{j,n}}{\lambda_{j,n+1}}\right)^{\epsilon}\right]\;.$
(26)

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

$\displaystyle{\boldsymbol{\lambda}_{n+1}}=\mathop{\mathrm{argmin}}_{{\boldsymbol{\lambda}}\in\mathcal{S}_{J}^{+}}\frac{1}{\epsilon}\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\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_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\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\ldots J$, we have by convexity of the function $u\mapsto u^{1+\epsilon}$
that
$\sum_{j=1}^{J}\left(\bar{\beta}_{j}/\lambda_{j}\right)^{1+\epsilon}\lambda_{j}\geqslant(\sum_{j=1}^{J}{\bar{\beta}_{j}})^{1+\epsilon}$
and that this quantity is minimal when $\lambda_{j}=\bar{\beta}_{j}$ for
$j=1\ldots 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_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\left[1-\left(\frac{\lambda_{j,n}}{\lambda_{j,n+1}}\right)^{\epsilon}\right]\leqslant
0$

which in turn yields (18) [since combined with (26) it implies (25) which
itself implies (18) as seen in the case $\eta_{n}=1$]. This settles the case
$\eta_{n}\in(0,1)$.

∎

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

${\boldsymbol{\lambda}_{n+1}}=\mathop{\mathrm{argmin}}_{{\boldsymbol{\lambda}}\in\mathcal{S}_{J}^{+}}h_{n}({\boldsymbol{\lambda}})$

where, setting $\epsilon=\eta_{n}^{-1}-1$, we have defined for all
${\boldsymbol{\lambda}}\in\mathcal{S}_{J}^{+}$,

$\displaystyle h_{n}({\boldsymbol{\lambda}})$
$\displaystyle=\begin{cases}\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\log\left(\frac{\lambda_{j}}{\lambda_{j,n}}\right),&\mbox{if
}\eta_{n}=1\;,\\\
\frac{1}{\epsilon}\sum_{j=1}^{J}\lambda_{j,n}\left[\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\nu(\mathrm{d}y)+\kappa_{n}\right]\left[1-\left(\frac{\lambda_{j,n}}{\lambda_{j}}\right)^{\epsilon}\right],&\mbox{if
}\eta_{n}\in(0,1)\;.\end{cases}$

More specifically, $h_{n}({\boldsymbol{\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 (23) in 3, we are
interested in deriving update formulas for the sequence
$(\Theta_{n})_{n\geqslant 1}$ satisfying (19), which we will then pair up with
(23) in order to apply 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\geqslant 1}$

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

####  4.2.1 A maximisation approach

As done in 1, an idea is to consider an update for $(\Theta_{n})_{n\geqslant
1}$ of the form:

$\displaystyle\Theta_{n+1}=\mathop{\mathrm{argmax}}_{\Theta\in{\mathsf{T}}^{J}}g_{n}(\Theta)\;,\quad
n\geqslant 1\;,$

where the function $g_{n}$ is constructed as a lower bound on
$\Theta\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$ of the function
$\Theta\mapsto\int_{\mathsf{Y}}\sum_{j=1}^{J}\lambda_{j,n}\varphi_{j,n}^{(\alpha)}(y)\log\left({k(\theta_{j},y)}/{k(\theta_{j,n},y)}\right)\nu(\mathrm{d}y)$
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})\leqslant 0$ implies (19). This leads us to 4
below.

######  Corollary 4 (Generalised maximisation approach).

Assume (A1). Let $\alpha\in[0,1)$, let $(\eta_{n})_{n\geqslant 1}$ be valued
in $(0,1]$ and let $(\kappa_{n})_{n\geqslant 1}$ be such that
$(\alpha-1)\kappa_{n}\geqslant 0$ at all times $n\geqslant 1$. Furthermore,
let $(b_{j,n})_{n\geqslant 1}$ be a non-negative sequence for all $j=1\ldots
J$. Starting from an initial parameter set
$({\boldsymbol{\lambda}_{1}},\Theta_{1})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$,
let $({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ be defined
iteratively for all $n\geqslant 1$ in such a way that (23) holds and

$\displaystyle\theta_{j,n+1}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{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(\mathrm{d}y)\;,\quad
j=1\ldots J\;{{,}}$ (27)

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

The proof of this result is deferred to Section B.1. Observe that, under the
assumptions of 4, Algorithm 1 leads to a systematic decrease in
$\Psi_{\alpha}$ at each step. This result effectively generalises the
monotonicity property from 1 to the case of mixture models as we can deduce
simple iterative schemes satisfying (27) for well-chosen kernel densities $k$.
To illustrate this, we provide below the extension of 1 to Gaussian Mixture
Models (and 1 can be extended in the same way, see Section B.2).

At iteration $n$,

For all $j=1\ldots J$, set  $\displaystyle\lambda_{j,n+1}$
$\displaystyle=\frac{\lambda_{j,n}\left[\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}{\sum_{\ell=1}^{J}\lambda_{\ell,n}\left[\int_{\mathsf{Y}}\varphi_{\ell,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}$
$\displaystyle\theta_{j,n+1}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{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(\mathrm{d}y)\;.$

Algorithm 1 Maximisation approach algorithm for mixture models

######  Example 3 (Maximisation approach for Gaussian Mixture Models).

We consider the case of $d$-dimensional Gaussian mixture densities, in which
case $\mathsf{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{\mathsf{T}}$ denotes the mean and covariance
matrix of the $j$-th Gaussian component density. We let the parameter set
${\mathsf{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 $\mathsf{Y}$ satisfies (8).

Then (A1) holds. Moreover, starting from any
$({\boldsymbol{\lambda}_{1}},\Theta_{1})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$
so that $\Psi_{\alpha}(\mu_{1}k)<\infty$ and using the notation
$\theta_{j,n}=(m_{j,n},\Sigma_{j,n})$ for all $j=1\ldots J$ and all
$n\geqslant 1$, it holds that: for all $j=1\ldots J$, all $n\geqslant 1$ and
all $\gamma_{j,n}\in(0,1]$, there exists $b_{j,n}\geqslant 0$ such that the
argmax problem (27) has a unique solution defined by:

$\displaystyle\quad m_{j,n+1}$
$\displaystyle=\gamma_{j,n}\widehat{m}_{j,n}+(1-\gamma_{j,n})m_{j,n}$
$\displaystyle\Sigma_{j,n+1}$
$\displaystyle=\gamma_{j,n}\widehat{\Sigma}_{j,n}+(1-\gamma_{j,n})\Sigma_{j,n}+\gamma_{j,n}(1-\gamma_{j,n})\left(\widehat{m}_{j,n}-m_{j,n}\right)\left(\widehat{m}_{j,n}-m_{j,n}\right)^{T}\;,$

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

$\displaystyle\widehat{m}_{j,n}=\int_{\mathsf{Y}}y\,\check{\varphi}_{j,n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)\;,$
$\displaystyle\widehat{\Sigma}_{j,n}=\int_{\mathsf{Y}}yy^{T}\,\check{\varphi}_{j,n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)-\widehat{m}_{j,n}\widehat{m}_{j,n}^{T}\;.$

The detailed derivation of the mean and covariance updates above is deferred
to Section 6.2. Notice that the interpretation of those updates is akin to the
one already made in 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 4, these updates are
compatible with the mixture weights updates (23), resulting in a systematic
decrease of $(\Psi_{\alpha}(\mu_{n}k))_{n\geqslant 1}$.

In line with 2, we next present another possible update formula for
$({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$.

####  4.2.2 A gradient-based approach

In the spirit of 2, we now resort to Gradient Descent steps to satisfy (18).

######  Corollary 5 (Generalised gradient-based approach).

Assume (A1). Let ${\mathsf{T}}\subseteq\mathbb{R}^{d}$ be a convex set, let
$\alpha\in[0,1)$, let $(\eta_{n})_{n\geqslant 1}$ be valued in $(0,1]$ and let
$(\kappa_{n})_{n\geqslant 1}$ be such that $(\alpha-1)\kappa_{n}\geqslant 0$
at all times $n$. Furthermore, for all $j=1\ldots J$, let
$(\gamma_{j,n})_{n\geqslant 1}$ be valued in $(0,1]$. Starting from an initial
parameter set
$({\boldsymbol{\lambda}_{1}},\Theta_{1})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$,
let $({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ be defined
iteratively for all $n\geqslant 1$ in such a way that (23) holds and

$\displaystyle\theta_{j,n+1}$
$\displaystyle=\theta_{j,n}-\frac{\gamma_{j,n}}{\beta_{j,n}}\nabla
g_{j,n}(\theta_{j,n})\;,\quad j=1\ldots J\;,$ (28)

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

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

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

At iteration $n$,

For all $j=1\ldots J$, set  $\displaystyle\lambda_{j,n+1}$
$\displaystyle=\frac{\lambda_{j,n}\left[\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}{\sum_{\ell=1}^{J}\lambda_{\ell,n}\left[\int_{\mathsf{Y}}\varphi_{\ell,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}$
$\displaystyle\theta_{j,n+1}$
$\displaystyle=\theta_{j,n}-\frac{\gamma_{j,n}}{\beta_{j,n}}\nabla
g_{j,n}(\theta_{j,n})\;.$

Algorithm 2 Gradient-based approach algorithm for mixture models

The proof of this result is deferred to Section B.3. Notice that, under the
assumptions of 5, Algorithm 2 ensures a systematic decrease in $\Psi_{\alpha}$
at each step and the framework of 5 thus extends 2 to the case of mixture
models. Much like what we already did for 2, an idea is then to identify how
our updates relate to Gradient Descent-based techniques for optimising
$\Theta$.

Under common differentiability assumptions, we can write: for all $n\geqslant
1$, all $j=1\ldots J$ and all $\theta^{\prime}\in{\mathsf{T}}$

$\nabla
g_{j,n}(\theta^{\prime})=\int_{\mathsf{Y}}\frac{\varphi_{j,n}^{(\alpha)}(y)}{\alpha-1}\left.\frac{\partial\log
k(\theta,y)}{\partial\theta}\right|_{(\theta,y)=(\theta^{\prime},y)}\nu(\mathrm{d}y)\;,$

so that (28) becomes

$\displaystyle\theta_{j,n+1}$
$\displaystyle=\theta_{j,n}-\frac{\gamma_{j,n}}{\beta_{j,n}}\int_{\mathsf{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(\mathrm{d}y)\;,\quad
j=1\ldots J\;.$ (29)

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 2 generalises to the case of Gaussian
Mixture Models.

Given
$({\boldsymbol{\lambda}_{n}},\Theta_{n})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$,
an index $j$ in $1\ldots J$ and letting $p=p(\cdot,\mathscr{D})$, performing
one Gradient Descent iteration w.r.t. $\theta_{j,n}$ for $\alpha$-divergence
and Rényi’s $\alpha$-divergence minimisation indeed respectively amounts to
updating $\theta_{j,n+1}$ as follows:

$\displaystyle\theta_{j,n+1}$
$\displaystyle=\theta_{j,n}-r_{j,n}\lambda_{j,n}\int_{\mathsf{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(\mathrm{d}y)\;,$
$\displaystyle\theta_{j,n+1}$
$\displaystyle=\theta_{j,n}-r_{j,n}\frac{\lambda_{j,n}}{\sum_{j=1}^{J}\lambda_{j,n}\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)}\int_{\mathsf{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(\mathrm{d}y)\;$

where $r_{j,n}>0$ is the learning rate (we refer to Section B.4 for details
regarding how these updates are derived). One can then see that there is a
similarity between (29) and the Gradient Descent updates above.

To fully comprehend the connection between our gradient-based approach and
Gradient Descent steps, we now present below the generalisation of 2 to
Gaussian Mixture Models. As we shall see, the smoothness assumption on
$g_{j,n}$ will be satisfied in that example and we will also 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 ${\mathsf{T}}\subseteq\mathbb{R}^{d}$
is a convex subset and that $\Sigma_{j}=\sigma_{j}^{2}\boldsymbol{I_{d}}$,
where $\sigma_{j}^{2}>0$ and $\boldsymbol{I_{d}}$ denotes the $d$-dimensional
identity matrix. Finally, we assume that the non-negative function $p$ defined
on $\mathsf{Y}$ satisfies (13).

Then, (A1) holds. Furthermore, setting
$\beta_{j,n}=\sigma_{j}^{-2}(1-\alpha)^{-1}\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)$,
the function $g_{j,n}$ is $\beta_{j,n}$-smooth for all $j=1\ldots J$ and all
$n\geqslant 1$. As a result, using the notation $\theta_{j,n}=m_{j,n}$, the
means of the Gaussian densities $(\mathcal{N}(y;m_{j},\Sigma_{j}))_{1\leqslant
j\leqslant J}$ can be optimised using the update formula at iteration
$n\geqslant 1$:

$\displaystyle m_{j,n+1}$
$\displaystyle=\gamma_{j,n}\hat{m}_{j,n}+\left(1-\gamma_{j,n}\right)m_{j,n}\;,\quad
j=1\ldots J\;,$

where $\check{\varphi}_{j,n}^{(\alpha)}$ is defined in (16),
$\gamma_{j,n}\in[0,1)$ and where we have set
$\hat{m}_{j,n}=\int_{\mathsf{Y}}y\leavevmode\nobreak\
\check{\varphi}_{j,n}^{(\alpha)}(y)\leavevmode\nobreak\ \nu(\mathrm{d}y)$ (we
refer to 8 for the detailed derivation).

In particular, letting the sequence $(\gamma^{\prime}_{j,n})_{n\geqslant 1}$
be valued in $(0,1]$ for all $j=1\ldots J$

$\gamma_{j,n}=\gamma_{j,n}^{\prime}\frac{\lambda_{j,n}\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)}{\int_{\mathsf{Y}}{\mu_{n}k(y)^{\alpha}p(y)^{1-\alpha}}\nu(\mathrm{d}y)}\in(0,1]$

(since
$\int_{\mathsf{Y}}\mu_{n}k(y)^{\alpha}p(y)^{1-\alpha}\nu(\mathrm{d}y)=\sum_{j=1}^{J}\lambda_{j,n}\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)$).
The resulting iterative algorithm is given by the following update at time
$n\geqslant 1$:

$\displaystyle m_{j,n+1}$
$\displaystyle=m_{j,n}+\gamma_{j,n}^{\prime}\frac{\int_{\mathsf{Y}}\lambda_{j,n}\varphi_{j,n}^{(\alpha)}(y)(y-m_{j,n})\nu(\mathrm{d}y)}{\int_{\mathsf{Y}}{\mu_{n}k(y)^{\alpha}p(y)^{1-\alpha}}\nu(\mathrm{d}y)}\;,\quad
j=1\ldots J\;.$ (30)

Interestingly, (30) 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}^{\prime}$. Hence, if we
were to solely rely on the Gradient Descent litterature, the convergence of
the iterative sequence $(\theta_{j,n})_{n\geqslant 1}$ defined by (30) would
require the sequence $({\boldsymbol{\lambda}_{n}})_{n\geqslant 1}$ to be
constant. This is in contrast with 5, which allows for a simultaneous
optimisation of ${\boldsymbol{\lambda}}$ and $\Theta$ according to (23) and
(30).

We can now add on the comments we made in Section 3 regarding the gradient-
based methodology built in 2 and that we have since extended to mixture models
in 5:

  * •

A core insight from 5, which is exemplified in 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 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 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
$({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ that meet the
conditions of 3 for a monotonic decrease in $\Psi_{\alpha}$ at each step. We
notably provided in Examples 3 and 4 simple updates formulas for Gaussian
Mixture Models, corresponding respectively to an application of the
maximisation approach from Section 4.2.1 and of the gradient-based approach
from Section 4.2.2. We move on to describing 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,
starting from the approach of [19].

### 5.1 Rényi Divergence Variational Inference [19]

In [19], 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\left\\{1\right\\}$,

$\displaystyle\mathcal{L}_{\alpha}(q,\mathscr{D}):=\frac{1}{1-\alpha}\log\left(\int_{\mathsf{Y}}q(y)^{\alpha}p(y,\mathscr{D})^{1-\alpha}\nu(\mathrm{d}y)\right)\;.$
(31)

Contrary to us, the work from [19] does not handle the case where $q$ belongs
to the mixture models variational family (15), that is
$q=\mu_{{\boldsymbol{\lambda}},\Theta}k$. Indeed, Gradient Ascent steps w.r.t.
${\boldsymbol{\lambda}}$ on (31) are not fit for handling the constrained
optimisation on the simplex that is required for mixture weights optimisation.

Yet, a parallel can be drawn in the GMM case between [19] and our approach by
observing that the gradient-based updates on the means in 4 each coincide with
a Gradient Ascent step on the objective function (31) 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).

This illustrates the fact that, in constrast with [19], our work provides a
theoretical framework which enables simultaneous optimisation of the mixture
weights ${\boldsymbol{\lambda}}$ and of the mixture components parameters
$\Theta$. In addition, note that 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.

We next move on to the Power Descent algorithm from [22].

### 5.2 The Power Descent algorithm [22]

In order to idenfity the connection between our work and the Power Descent
algorithm introduced in [22], let us briefly present the latter.

The Power Descent is a gradient-based algorithm that operates on probability
measures and performs $\alpha$-divergence minimisation for all
$\alpha\in\mathbb{R}\setminus\left\\{1\right\\}$. More precisely, equipping
${\mathsf{T}}$ with a $\sigma$-field $\mathcal{T}$ and denoting by
$\mathrm{M}_{1}({\mathsf{T}})$ the space of probability measures on
$({\mathsf{T}},\mathcal{T})$, the Power Descent optimises $\Psi_{\alpha}(\mu
k)$ with respect to $\mu\in\mathrm{M}_{1}({\mathsf{T}})$, where $\mu
k(y)=\int_{\mathsf{T}}\mu(\mathrm{d}\theta)k(\theta,y)$ for all
$\mu\in\mathrm{M}_{1}({\mathsf{T}})$ and all $y\in\mathsf{Y}$. Given an
initial probability measure $\mu_{1}\in\mathrm{M}_{1}({\mathsf{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\geqslant 1\;,$ (32)

where, for all $\mu\in\mathrm{M}_{1}({\mathsf{T}})$ and all
$\theta\in{\mathsf{T}}$,

$\displaystyle{b_{\mu,\alpha}}(\theta)$
$\displaystyle={\int_{\mathsf{Y}}}k(\theta,y)\frac{1}{\alpha-1}\left[\left(\dfrac{\mu
k(y)}{p(y)}\right)^{\alpha-1}-1\right]\nu(\mathrm{d}y)$
$\displaystyle\mathcal{I}_{\alpha}(\mu)(\mathrm{d}\theta)$
$\displaystyle=\dfrac{\mu(\mathrm{d}\theta)\leavevmode\nobreak\
[(\alpha-1)({b_{\mu,\alpha}}(\theta)+\kappa)+1]^{\frac{\eta}{1-\alpha}}}{\mu([(\alpha-1)({b_{\mu,\alpha}}+\kappa)+1]^{\frac{\eta}{1-\alpha}})}\;.$

Interestingly, a monotonicity result has been proved for the Power Descent
algorithm in [22]. In this result, that is obtained as a particular case of
[22, Theorem 1] with $\Gamma=[(\alpha-1)v+1]^{\eta/(1-\alpha)}$, the monotonic
decrease in $\Psi_{\alpha}$ of the scheme (32) holds for all
$\mu\in\mathrm{M}_{1}({\mathsf{T}})$, all
$\alpha\in\mathbb{R}\setminus\left\\{1\right\\}$, all $\eta\in(0,1]$ and all
$\kappa$ such that $(\alpha-1)\kappa\geqslant 0$. We provide below a slightly
more general version of this result, in which the monotonic decrease in
$\Psi_{\alpha}$ also holds for well-chosen values of $\eta$ that are larger
than $1$ when $\alpha<0$ (we refer to Section 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.

  1. (i)

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

  2. (ii)

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

  3. (iii)

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

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

  1. (i)

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

  2. (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 [22, Theorem 1] for the Power
Descent algorithm - that we generalised in 2 \- [22] then applied this
algorithm to mixture weights optimisation by letting the initial probability
measure $\mu_{1}\in\mathrm{M}_{1}({\mathsf{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{\mathsf{T}}$ and
${\boldsymbol{\lambda}_{1}}\in\mathcal{S}_{J}$. For that choice of $\mu_{1}$,
$\mu_{n}$ in (32) 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:

$\displaystyle\lambda_{j,n+1}$
$\displaystyle=\dfrac{\lambda_{j,n}\leavevmode\nobreak\
[(\alpha-1)({b_{\mu_{n},\alpha}}(\theta_{j})+\kappa)+1]^{\frac{\eta}{1-\alpha}}}{\sum_{\ell=1}^{J}\lambda_{\ell,n}[(\alpha-1)({b_{\mu_{n},\alpha}}(\theta_{\ell})+\kappa)+1]^{\frac{\eta}{1-\alpha}}}\;,\quad
j=1\ldots J.$ (33)

One may then observe that the update (33) corresponds to the update on the
mixture weights (23) we have identified in 3 for $\Theta_{n}=\Theta$,
$\eta_{n}=\eta/(1-\alpha)$ and $\kappa_{n}=\kappa$. Steaming from this link
between our approach and [22], we now make two important comments.

  * •

Benefits of our approach compared to [22]. The monotonicity result from [22,
Theorem 1], that is generalised in 2, requires the sequence
$(\Theta_{n})_{n\geqslant 1}$ to be constant when applied to mixture models
optimisation as done in (33). This restricts the variational family to mixture
models with fixed mixture components parameters.

To remedy this problem, [22] proposed a fully-adaptive algorithm that
alternates between an Exploitation step optimising the mixture weights
according to (33) and an Exploration step acting on the mixture components
parameters. However, they did not provide choices of Exploration steps that
optimise the mixture components parameters and therefore, there is no
theoretical guarantees for their complete Exploitation-Exploration algorithm.

A strong improvement of our approach is then that we provide updates for
$({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ that systematically
decrease $\Psi_{\alpha}$ at each step, with the particularity that our mixture
weights updates relate to the framework of [22]. In that sense, our approach
supplements the work done in [22] (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. [22] 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 (33) 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 (23) from 3
and provides a better understanding of the role of the parameter $\eta_{n}$
appearing in this update. This aspect will notably come in handy when
interpreting our numerical experiments in Section 8.

We wrap up Section 5 by investigating the connections between our framework
and the Mixture-Population Monte Carlo (M-PMC) algorithm from [23].

### 5.3 The M-PMC algorithm [23]

For any measurable positive function $p$ on $(\mathsf{Y},\mathcal{Y})$, the
M-PMC algorithm [23] aims at solving the optimisation problem

$\displaystyle\sup_{({\boldsymbol{\lambda}}\in\mathcal{S}_{J},\Theta\in{\mathsf{T}}^{J})}\int_{\mathsf{Y}}\log\left(\sum_{j=1}^{J}\lambda_{j}k(\theta_{j},y)\right)p(y)\nu(\mathrm{d}y)\;,$

or equivalently, using a Variational Inference formulation, at minimising the
inclusive Kullback-Leibler divergence
$D_{0}(\mu_{{\boldsymbol{\lambda}},\Theta}K||\mathbb{P})$ w.r.t. to
$({\boldsymbol{\lambda}},\Theta)$, where
$\mathbb{P}(A)=\int_{A}p\leavevmode\nobreak\
\mathrm{d}\nu/\int_{\mathsf{Y}}p\leavevmode\nobreak\ \mathrm{d}\nu$ for all
$A\in\mathcal{Y}$. This is done in [23, Section 2] by introducing the
following iterative update formulas: for all $j=1\ldots J$ and all $n\geqslant
1$

$\displaystyle\lambda_{j,n+1}$
$\displaystyle=\int_{\mathsf{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_{\mathsf{Y}}p(y)\nu(\mathrm{d}y)}\nu(\mathrm{d}y)$
$\displaystyle\theta_{j,n+1}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta_{j}\in{\mathsf{T}}}\int_{\mathsf{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(\mathrm{d}y)\;.$

Interestingly, [23] noticed that the two updates above 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

$\displaystyle\sup_{({\boldsymbol{\lambda}}\in\mathcal{S}_{J},\Theta\in{\mathsf{T}}^{J})}\sum_{m=1}^{M}\log\left(\sum_{j=1}^{J}\lambda_{j}k(\theta_{j},Y_{m})\right)\;.$

Observe now that those update formulas correspond to having considered the
particular case $\alpha=0$, $\eta_{n}=1$, $\kappa_{n}=0$, $a_{j,n}=1$ and
$b_{j,n}=0$ in Algorithm 1, meaning that the M-PMC algorithm is a special case
of our framework. As we shall see, the additionnal layer of flexibility we
have obtained in Algorithm 1 will 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.

Furthermore, we can interpret the results we have obtained thus far in the
light of the M-PMC algorithm: we have found a generalisation of an integrated
EM algorithm preserving the monotonicity property for a wide range of
hyperparameters. Here, a particularly striking fact is that the monotonicity
property holds for $\alpha\in[0,1)$, so that updates akin to an Integrated EM
procedure can be derived past the traditional case $\alpha=0$ corresponding to
integrated likelihood optimisation.

A table summarising the main improvements of our framework compared to the
existing litterature can be found in Table 1 and in the coming section, we
revisit the maximisation and the gradient-based approaches when the
variational family is based on the exponential family.

| Improvements of our framework
---|---
Rényi Gradient Descent | Simultaneous optimisation w.r.t. $({\boldsymbol{\lambda}},\Theta)$
[19] | (prev. mixture weights ${\boldsymbol{\lambda}}$ optimisation not considered)
| For GMMs : maximisation approach encompasses Rényi
| Gradient Descent and provides covariance matrices updates
Power Descent | Simultaneous optimisation w.r.t. $({\boldsymbol{\lambda}},\Theta)$
[22] | (prev. $(\Theta_{n})_{n\geqslant 1}$ constant)
M-PMC algorithm | Extension of an Integrated EM algorithm to:
[23] |  $\alpha\in[0,1)$, $\eta_{n}\in(0,1]$, $(\alpha-1)\kappa_{n}\geqslant 0$ and $b_{j,n}\geqslant 0$
| (prev. $\alpha=0$, $\eta_{n}=1$, $\kappa_{n}=0$ and $b_{j,n}=0$)
Table 1: Key improvements of our framework for mixture models optimisation
compared to related works (prev. stands for previously in the literature).

## 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.

######  Definition 1 (Exponential family).

Let $\nu$ be a $\sigma$-finite measure on $(\mathsf{Y},\mathcal{Y})$, $E$ be a
Euclidean space endowed with an inner product
$\left\langle\cdot,\cdot\right\rangle_{E}$ and its corresponding norm
$\left\|\cdot\right\|_{E}$, $h$ be a non-negative function defined on
$\mathsf{Y}$ and $S:\mathsf{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\left(\left\langle\zeta,S(y)\right\rangle_{E}-A(\zeta)\right)\;,\quad\;y\in\mathsf{Y}\;,\;\zeta\in
E_{0}$

where

$E_{0}:=\left\\{\zeta\in
E\;:\;\int_{\mathsf{Y}}h(y)\;\mathrm{e}^{\left\langle\zeta,S(y)\right\rangle_{E}}\;\nu(\mathrm{d}y)<\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)\;\mathrm{d}\nu=1$. In this
setting, $\nu$ is called the dominating measure and $S$ the natural statistic.
In what follows, $\mathrm{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
$\mathrm{Int}E_{0}$. Moreover, for all $\zeta\in\mathrm{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\nabla^{T}A$ of $A$ taken at $\zeta$, respectively. For later
reference, we write in particular

$\displaystyle\nabla
A(\zeta)=\int_{\mathsf{Y}}S(y)\,k^{(o)}(\zeta,y)\;\nu(\mathrm{d}y)\;$ (34)

and we introduce the following assumption.

  1. (B1)

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

    1. (a)

The function $h$ on $\mathsf{Y}$ is positive.

    2. (b)

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

Notice that the assumption (B1)(B1)(a) can be circumvented by changing the
dominating measure $\nu$ and therefore the assumptions (B1)(B1)(a) and
(B1)(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\mathrm{Int}E_{0}$. In addition, one often relies on a non-canonical
parameterisation (see for example 1) leading to the kernel density

$\displaystyle
k(\theta,y)=k^{(o)}(\upsilon(\theta),y)=h(y)\;\exp\left(\left\langle\upsilon(\theta),S(y)\right\rangle_{E}-A\circ\upsilon(\theta)\right)\;,\quad\theta\in{\mathsf{T}}\;,\;y\in\mathsf{Y}\;,$
(35)

where $\upsilon$ maps the parameter space ${\mathsf{T}}$ to a subset of
$\mathrm{Int}E_{0}$. The assumption (B1)(B1)(a) then ensures that the
assumption made on $k$ in (A1) is satisfied, which will prove to be useful in
the upcoming results.

### 6.2 The maximisation approach for the exponential family

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

$\theta_{n+1}=\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{Y}}\left[\ratio(y)+b_{n}k(\theta_{n},y)\right]\log\left(\frac{k(\theta,y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\;,\quad
n\geqslant 1\;.$

Furthermore, we have seen that in the more general case of mixture models
described in 4, the maximisation approach leads to the argmax problem (27).
This second argmax problem is similar to (7) in the sense that it updates
$\Theta_{n}=(\theta_{j,n})_{1\leqslant j\leqslant J}$ by solving $J$
component-wise argmax problems of the same form as (7) (with $\theta_{n+1}$, ,
$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\dots J$):

$\theta_{j,n+1}=\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{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(\mathrm{d}y)\;,\quad
j=1\ldots J\;,\;n\geqslant 1\;.$

Our goal is now to solve the argmax problems above for a member of the
exponential family. To do so, we will first state a theorem in which the
kernel $k$ is in its canonical form, that is $v$ is the identity mapping in
(35) and $k(\theta,\cdot)=k^{(o)}(\zeta,\cdot)$.

Other parameterisations will then be built on this case and as a result we
will deduce corollaries that include non-canonical parameterisations and that
are applicable to the general case of mixture models. Note that all the proofs
of the results stated in Section 6.2 are deferred to Section 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\left\|S\right\|_{E}\,\check{\varphi}\;\mathrm{d}\nu<\infty\;.$ (36)

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

$\mathop{\mathrm{argmax}}_{\zeta\in\mathcal{O}}\int_{\mathsf{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(\mathrm{d}y)$
(37)

admits at least a solution if and only if

$\mathbf{s}^{*}:=\frac{1}{1+b}\int
S\,\check{\varphi}\;\mathrm{d}\nu+\frac{b}{1+b}\nabla A(\zeta_{0})$ (38)

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

$\nabla A(\zeta^{*})=\mathbf{s}^{*}\;.$ (39)

4 shows the equivalence between the argmax problem (37) and the equation (39)
under the canonical parameterisation. From there, using the expression of
$\nabla A$ given in (34) and considering the non-canonical exponential family
probability density (35), we can interpret the equation (39) as an equality
between two means of the statistic $S$ computed under two different
distributions. More specifically, setting $\zeta^{*}=\upsilon(\theta^{*})$ and
$\zeta_{0}=\upsilon(\theta_{0})$, we obtain

$\int_{\mathsf{Y}}S(y)k(\theta^{*},y)\nu(\mathrm{d}y)=\int_{\mathsf{Y}}S(y)\psi(y)\nu(\mathrm{d}y)$
(40)

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\mathsf{Y}\;.$

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

$\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{Y}}\left[\check{\varphi}(y)+bk(\theta_{0},y)\right]\log\left(\frac{k(\theta,y)}{k(\theta_{0},y)}\right)\nu(\mathrm{d}y)\;.$
(41)

######  Corollary 6 (Gaussian density).

Let $d\geqslant 1$ and $\mathsf{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{\mathsf{T}}=\mathbb{R}^{d}\times\mathcal{M}_{>0}(d)$.
Let $\check{\varphi}$ be a probability density function on $\mathsf{Y}$ such
that (36) holds, that is in the Gaussian case

$\int_{\mathsf{Y}}\|y\|^{2}\,\check{\varphi}(y)\;\mathrm{d}y<\infty\;.$

Let $\theta_{0}=(m_{0},\Sigma_{0})\in{\mathsf{T}}$ and $b\geqslant 0$. If
$b=0$, assume moreover that $\check{\varphi}$ is non-degenerate in the sense
that its covariance matrix is positive definite. Then (41) has a unique
solution $\theta^{*}=(m^{*},\Sigma^{*})\in{\mathsf{T}}$ defined by

$\displaystyle m^{*}=\mathbb{E}[Y]$
$\displaystyle\Sigma^{*}=\mathbb{C}\mathrm{ov}(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})\;.$

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

######  Corollary 7.

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

$0<\int_{\mathsf{Y}}\left(1+\|S(y)\|_{E}^{1/(1-\alpha)}\right)\;p(y)\;\nu(\mathrm{d}y)<\infty\;$
(42)

holds. Let $J\in\mathbb{N}$, $n\geqslant 1$ and
$({\boldsymbol{\lambda}_{n}},\Theta_{n})\in\mathcal{S}_{J}^{+}\times{\mathsf{T}}^{J}$.
Then, for all $j=1\ldots J$, (36) holds with
$\check{\varphi}=\check{\varphi}_{j,n}^{(\alpha)}$ as defined in (16) and
setting
$\gamma_{j,n}=\int\varphi_{j,n}^{(\alpha)}\mathrm{d}\nu/(\int\varphi_{j,n}^{(\alpha)}\mathrm{d}\nu+b_{j,n})$,
we have that

$\displaystyle\mathop{\mathrm{argmax}}_{\theta\in{\mathsf{T}}}\int_{\mathsf{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(\mathrm{d}y)$
(43)

admits at least a solution if and only if

$\displaystyle\mathbf{s}_{j}^{*}:=\gamma_{j,n}\int S\leavevmode\nobreak\
\check{\varphi}_{j,n}^{(\alpha)}\leavevmode\nobreak\
\mathrm{d}\nu+(1-\gamma_{j,n})\nabla A\circ\upsilon(\theta_{j,n})$ (44)

belongs to the image set $\nabla A\circ\upsilon({\mathsf{T}})$, in which case
(43) has a unique solution $\theta_{j}^{*}$ defined by

$\displaystyle\nabla A\circ\upsilon(\theta_{j}^{*})=\mathbf{s}_{j}^{*}$

and we can set $\theta_{j,n+1}=\theta^{*}_{j}$ for all $j=1\ldots J$ in (27)
of 4.

Let us now comment on 7.

  * •

The updates from Examples 1 and 3 follow from 7, paired up with 6 and with the
fact that the condition (42) simplifies to (8) (since in the Gaussian case we
have set $S(y)=(y,-yy^{T}/2)$ for all $y\in\mathsf{Y}$).

  * •

More generally, we can always solve the argmax problem (43), that is the
argmax problem of 4 with $k$ belonging to the exponential family, for
$b_{j,n}>0$ large enough.

Indeed, under the assumptions of 7, $\upsilon({\mathsf{T}})$ is an open subset
of $E_{0}$ and $\nabla A$ is a $\mathcal{C}^{\infty}$-diffeomorphism on
$\mathrm{Int}E_{0}$. As a result, $\nabla A\circ\upsilon({\mathsf{T}})$ is an
open subset of $E$. Since $\mathbf{s}_{j}^{*}$ in (44) tends to $\nabla
A\circ\upsilon(\theta_{j,n})$ as $b_{j,n}\to\infty$, we get that
$\mathbf{s}_{j}^{*}$ always ends up belonging to $\nabla
A\circ\upsilon({\mathsf{T}})$ for $b_{j,n}$ large enough. Furthermore, recall
that the regularisation parameter $b_{j,n}>0$ in 4 can be chosen at each step
$n$ of the iterative algorithm. Hence, under the assumptions of 7, we can
always solve the argmax problem of 4 for $b_{j,n}>0$ large enough.

  * •

Thanks to its argmax form, the argmax problem (43) is immune to change of
variables. The maximisation approach is thus solved using the canonical
parameter $\zeta$ in 4 and then transported to the parameter $\theta$ via the
one-to-one mapping $v$ as written in 7.

In contrast and as we shall see next in Section 6.3, the mapping $v$ will have
a strong impact on the gradient-based approach.

### 6.3 The gradient-based approach for an exponential family distribution

Remember that the key aspect to construct the sequence
$(\theta_{n})_{n\geqslant 1}$ according to the gradient-based approach from 2
lies in (i) computing the gradient of the function $g_{n}$ defined in (11) by:
for all $\theta\in{\mathsf{T}}$,

$\displaystyle
g_{n}(\theta)=\int_{\mathsf{Y}}\frac{{{\ratio(y)}}}{\alpha-1}\log\left(\frac{k(\theta,y)}{k(\theta_{n},y)}\right)\nu(\mathrm{d}y)\;$

and (ii) verifying that the function $g_{n}$ satisfies a
$\beta_{n}$-smoothness condition, in order to perform updates of the form

$\theta_{n+1}=\theta_{n}-\frac{\gamma_{n}}{\beta_{n}}\nabla
g_{n}(\theta_{n})\;,\quad n\geqslant 1$

with $\gamma_{n}\in(0,1]$. As seen in 5, this procedure is extended to mixture
models in a way that is similar to 2: it involves the gradient of $J$
functions $g_{j,n}$ that resemble $g_{n}$ ($g_{n}$, and $\theta_{n}$ are
replaced by $g_{j,n}$, $\varphi_{j,n}^{(\alpha)}$ and $\theta_{j,n}$ in the
definition above) and where each function $g_{j,n}$ is assumed to be
$\beta_{j,n}$-smooth.

We are now interested in solving the gradient-based approach for a member of
the exponential family. Like in the previous section, we will start with a
result handling the case where the kernel $k$ is in its canonical form (that
is, $v$ is the identity mapping in (35) and
$k(\theta,\cdot)=k^{(o)}(\zeta,\cdot)$) and we will defer the proofs of
Section 6.3 to Section C.3.

######  Proposition 3.

Let $k^{(o)}$ satisfy (B1). Let $\check{\varphi}$ be a probability density
function with respect to $\nu$ satisfying (36), that is

$\int\left\|S\right\|_{E}\,\check{\varphi}\;\mathrm{d}\nu<\infty\;.$

Let $\zeta_{0}\in\mathrm{Int}E_{0}$ and set

$g^{(o)}(\zeta)=-\int_{\mathsf{Y}}\check{\varphi}(y)\;\log\left(\frac{k^{(o)}(\zeta,y)}{k^{(o)}(\zeta_{0},y)}\right)\nu(\mathrm{d}y)\;,\quad\zeta\in\mathrm{Int}E_{0}\;.$
(45)

Then, for any convex subset $C_{0}\subseteq\mathrm{Int}E_{0}$, the function
$g^{(o)}$ is $\beta_{0}$-smooth over $C_{0}$ with

$\displaystyle\nabla g^{(o)}(\zeta)=\nabla A(\zeta)-\int
S\,\check{\varphi}\;\mathrm{d}\nu\quad\text{and}\quad\beta_{0}\geqslant\sup_{\zeta\in
C_{0}}\left\|\nabla\nabla^{T}A(\zeta)\right\|_{\mathrm{op}}\;,$ (46)

where $\left\|\cdot\right\|_{\mathrm{op}}$ denotes the operator norm.

Observe that, letting $C_{0}$ be the line segment
$\\{\zeta_{0}-\frac{t}{\beta_{0}}\left(\nabla A(\zeta_{0})-\int
S\,\check{\varphi}\;\mathrm{d}\nu\right)\;:\;{t\in[0,1]}\\}$ in 3, we have
that $C_{0}$ stays in $\mathrm{Int}E_{0}$ for $\beta_{0}$ large enough. Since
the Hessian of $A$ is continuous on $\mathrm{Int}E_{0}$, the inequality in
(46) can also be satisfied. Hence, for all $\gamma\in(0,1]$, the following
gradient step

$\displaystyle\zeta=\zeta_{0}-\frac{\gamma}{\beta_{0}}\left(\nabla
A(\zeta_{0})-\int_{\mathsf{Y}}S\,\check{\varphi}\;\mathrm{d}\nu\right)\;$ (47)

can always be performed and as a consequence we can always apply the gradient-
based approach when the kernel $k$ is in its canonical form.

As a result, when considering non-canonical exponential family probability
densities of the form (35) and under some reasonable assumptions on the
mapping $\upsilon$ (for instance, if $\upsilon$ is a $\mathcal{C}^{2}$
diffeomorphism), one is left with choosing between two alternative ways of
exploiting 3:

  1. (i)

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

$\displaystyle\theta=\upsilon^{-1}\left(\upsilon(\theta_{0})-\frac{\gamma}{\beta_{0}}\left(\nabla
A\circ\upsilon(\theta_{0})-\int_{\mathsf{Y}}S\,\check{\varphi}\;\mathrm{d}\nu\right)\right)\;.$

  2. (ii)

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

$g(\theta)=-\int\check{\varphi}(y)\;\log\left(\frac{k(\theta,y)}{k(\theta_{0},y)}\right)\mathrm{d}y\;,\quad\theta\in{\mathsf{T}}\;.$

In that case, since $g=g^{(o)}\circ\upsilon$ where $g^{(o)}$ is again as in
(45) with $\zeta_{0}=\upsilon(\theta_{0})$, this leads to an update of the
form

$\displaystyle\theta=\theta_{0}-\frac{\gamma}{\beta}\nabla\upsilon(\theta_{0})\cdot\left(\nabla
A\circ\upsilon(\theta_{0})-\int_{\mathsf{Y}}S\,\check{\varphi}\;\mathrm{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\geqslant 1$ and $\mathsf{Y}=\mathbb{R}^{d}$. We consider the
exponential family of a $d$-dimensional Gaussian density with
$k(\theta,y)=\mathcal{N}(y;\theta,\Sigma)$, where
$\theta\in{\mathsf{T}}:=\mathbb{R}^{d}$ and $\Sigma$ is a (known) covariance
matrix in $\mathcal{M}_{>0}(d)$. Let $\check{\varphi}$ be a probability
density function with respect to $\nu$ such that (36) holds, that is in this
particular case

$\int\|y\|\,\check{\varphi}(y)\;\mathrm{d}y<\infty\;.$

Set $\upsilon(\theta)=\Sigma^{-1}\theta$ and define the canonical kernel
$k^{(o)}$ with
$h(y)=(2\pi)^{-d/2}\left|\Sigma\right|^{-1/2}\,\mathrm{e}^{-y^{T}\Sigma^{-1}y/2}$,
$S(y)=y$ and $A(\zeta)=\zeta^{T}\Sigma\zeta/2$. Then (35) holds and Methods
(i) and (ii) lead to the the following gradient-based updates, respectively,

$\displaystyle\theta=\theta_{0}-\frac{\gamma}{\beta_{0}}\Sigma\left(\theta_{0}-\int
y\leavevmode\nobreak\ \check{\varphi}(y)\;\mathrm{d}y\right)\;,$ (48)
$\displaystyle\theta=\theta_{0}-\frac{\gamma}{\beta}\Sigma^{-1}\,\left(\theta_{0}-\int
y\check{\varphi}(y)\;\mathrm{d}y\right)\;,$ (49)

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.

Crucially, 8 illustrates how, even in a simple framework, the gradient-based
update formulas strongly depend on the parameterisation chosen. Indeed, they
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 8 with
$\check{\varphi}=\check{\varphi}_{n}^{(\alpha)}$ (and more generally with
$\check{\varphi}=\check{\varphi}_{j,n}^{(\alpha)}$ for Gaussian Mixture
Models). As the condition (36) is implied by (13) (using 8 in Section C.2 and
that here $S(y)=y$), this enables us to deduce the update formulas in Examples
2 and 4. Note that we did not introduce a convex subset $C_{0}$ in 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 2 and 4.

To conclude this part, let us put into perspective 6, in which we considered a
Gaussian density with varying mean vector and covariance matrix, and 8, where
the covariance matrix is assumed to be known.

While the simplified framework of 8 allowed us to keep 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
exactly the same setting as in 6. Observe indeed that for
$\theta=(m,\Sigma)\in{\mathsf{T}}=\mathbb{R}^{d}\times\mathcal{M}_{>0}(d)$ and
under the condition (36), the gradient of $g$ in that specific case is given
by

$\displaystyle\nabla g(\theta)$ $\displaystyle=\left(\Sigma^{-1}\left(m-\int
y\check{\varphi}(y)\;\mathrm{d}y\right),\frac{1}{2}\Sigma^{-1}\left(\int
yy^{T}\check{\varphi}(y)\;\mathrm{d}y-\Sigma+mm^{T}\right)\Sigma^{-1}\right)\;.$

Unfortunately, the smoothness index $\beta$ is now no longer constant, which
makes it much more involved to choose a convenient convex subset $C_{0}$ and
to bound the smoothness index over it.

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

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 $\Sigma$ and $a$ degrees of freedom which can be
defined as the continuous mixture

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

where $\check{\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. Note that
all the proofs from Section 6.4 are deferred to Section C.4.

####  6.4.1 The (E)LM family: definition and notation

Let $F$ be an Euclidean space endowed with the inner product
$\left\langle\cdot,\cdot\right\rangle_{F}$ and norm
$\left\|\cdot\right\|_{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 then provide below
the definition of the (E)LM family.

######  Definition 2 ((E)LM family).

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

  1. (i)

for all $\tau,\tau^{\prime}$ in $\mathrm{M}$, the distributions $\tau$ and
$\tau^{\prime}$ are equivalent;

  2. (ii)

for all $\tau\in\mathrm{M}$ and $\tau$-almost all $\ell\in\mathsf{Z}$, the
image set $\ell\circ\Upsilon(\tilde{{\mathsf{T}}})$ is included in $E_{0}$.

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

$\displaystyle
k{{{}^{(1)}}}(\theta,y)=\int_{\mathsf{Z}}k^{(o)}(\ell\circ\Upsilon(\vartheta),y)\;\tau(\mathrm{d}\ell)\;,\qquad
y\in\mathsf{Y}\;.$ (50)

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

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

The mapping $\Upsilon$ in 2 allows us to introduce various parameterisations.
However it is important to note that for
$\theta=(\vartheta,\tau)\in{\mathsf{T}}$, the probability density function
$k^{(1)}(\theta,\cdot)$ only depends on
$\xi=\Upsilon(\vartheta)\in\Upsilon(\tilde{{\mathsf{T}}})$ and
$\tau\in\mathrm{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 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 (27) as we did for the exponential family in Section 6.2. Instead, we
come back to a monotonic decrease condition of the form

$\int_{\mathsf{Y}}\check{\varphi}(y)\log\left(\frac{k^{(1)}(\theta,y)}{k^{(1)}(\theta_{0},y)}\right)\nu(\mathrm{d}y)\geqslant
0\;,$ (51)

for an LM family $k^{(1)}$ as in 2. For selected choices of $\check{\varphi}$
and of $(\theta,\theta_{0})$, this condition can indeed be linked to the
Conditions (5) and (19) appearing in 1 and 2 respectively. This is detailed in
the remark below.

######  Remark 1.

Let $\alpha\in[0,1)$. Taking $\check{\varphi}=\check{\varphi}_{n}^{(\alpha)}$
defined by (6) and $(\theta,\theta_{0})$ replaced by
$(\theta_{n+1},\theta_{n})$, Condition (51) becomes Condition (5) of 1.
Applying Condition (51) for $j=1,\dots,J$ with
$\check{\varphi}=\check{\varphi}_{j,n}^{(\alpha)}$ defined by (16) and
$(\theta,\theta_{0})$ replaced by $(\theta_{j,n+1},\theta_{j,n})$ yields
Condition (19) of 2.

Now let
$\mathrm{D}_{\tau^{\prime},\tau}=\frac{\mathrm{d}\tau^{\prime}}{\mathrm{d}\tau}$
denotes the Radon-Nikodym derivative of $\tau^{\prime}$ with respect to
$\tau$. We then have the following proposition.

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

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

$\displaystyle\int_{\mathsf{Y}}\left(\int_{\mathsf{Z}}\varphi(\ell,y)\;\log\left(\mathrm{D}_{\tau,\tau_{0}}(\ell)\right)\tau_{0}(\mathrm{d}\ell)\right)\nu(\mathrm{d}y)\geqslant
0\;,$ (52)
$\displaystyle\int_{\mathsf{Y}}\left(\int_{\mathsf{Z}}\varphi(\ell,y)\log\left(\frac{k^{(o)}(\ell\circ\Upsilon(\vartheta),y)}{k^{(o)}(\ell\circ\Upsilon(\vartheta_{0}),y)}\right)\tau_{0}(\mathrm{d}\ell)\right)\nu(\mathrm{d}y)\geqslant
0$ (53)

where we set

$\displaystyle\varphi(\ell,y)=\check{\varphi}(y)\,\frac{k^{(o)}(\ell\circ\Upsilon(\vartheta_{0}),y)}{k^{(1)}(\theta_{0},y)}\;.$
(54)

Then Condition (51) holds with $\theta=(\vartheta,\tau)$ and
$\theta_{0}=(\vartheta_{0},\tau_{0})$ .

One possible way to obtain (52) is to observe the following.

######  Remark 2.

Since the left-hand side of (52) is zero for $\tau=\tau_{0}$, Condition (52)
is fulfilled by setting

$\displaystyle\tau=\mathop{\mathrm{argmax}}_{\tau^{\prime}\in\mathrm{M}}\int_{\mathsf{Y}}\left(\int_{\mathsf{Z}}\varphi(\ell,y)\;\log\left(\mathrm{D}_{\tau^{\prime},\tau_{0}}(\ell)\right)\tau_{0}(\mathrm{d}\ell)\right)\nu(\mathrm{d}y)\;,$
(55)

assuming that this argmax is well-defined. This will be done in order to get
the update formula (72) of 5, see the proof of 10 in Section C.4.

As for the updating of $\vartheta$, we can again adopt a maximisation approach
to derive a convenient $\vartheta$, as in the following result, where we
directly update $\xi$, since in 2, $\vartheta=\Upsilon(\xi)$ with a known
mapping $\Upsilon$.

######  Corollary 9.

Consider an LM family as in 2 and
$\theta_{0}=(\vartheta_{0},\tau_{0})\in\tilde{{\mathsf{T}}}\times\mathrm{M}$.
Let $\check{\varphi}:\mathsf{Y}\to\mathbb{R}_{+}$ be a probability density
function with respect to $\nu$ and define
$\varphi:\mathsf{Z}\times\mathsf{Y}\to\mathbb{R}_{+}$ by (54). Let
$\vartheta\in\tilde{{\mathsf{T}}}$ such that $\Upsilon(\vartheta)$ is a
solution of the argmax problem

$\displaystyle\mathop{\mathrm{argmax}}_{\xi\in\Upsilon(\tilde{{\mathsf{T}}})}\int_{\mathsf{Y}}\left(\int_{\mathsf{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}(\mathrm{d}\ell)\right)\nu(\mathrm{d}y)\;,$
(56)

where $\xi_{0}=\Upsilon(\vartheta_{0})$ and $b_{0}$ is any non-negative
function defined on $\mathsf{Z}$ such that $\int
b_{0}\;\mathrm{d}\tau_{0}<\infty$. Then Condition (53) holds.

Interestingly, the maximisations (55) and (56) can be achieved for continuous
mixtures such as the one used to define Student’s $t$ distribution but they
require a thorough analysis. This is the goal of the next section, where we
formulate and solve a maximisation approach for a kernel $k^{(1)}$ that
belongs to the ELM family.

####  6.4.3 The maximisation approach for the ELM family

In the next result, an argmax of the form (56) is solved 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 (hence satisfies (B1))
and by relying on additional assumptions to be introduced alongside with some
helpful notation.

Our first assumption is the following.

  1. (C1)

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

Assumption (C1) guaranties that, for all $\vartheta\in\tilde{{\mathsf{T}}}$,
$\tau\in\mathrm{M}$ and $\tau$-almost all $\ell$ in $\mathsf{Z}$,
$\ell\circ\Upsilon(\vartheta)$ belongs to $\mathrm{Int}E_{0}$, the interior
set of $E_{0}$ (see 9 in Section 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
$\left\|\cdot\right\|_{\mathrm{op}}$ denote operator norms, for instance, for
any $\ell\in\mathsf{Z}$,

$\left\|\ell\right\|_{\mathrm{op}}=\sup\left\\{\left\|\ell(x)\right\|_{F}\;:\;x\in
E\,,\,\left\|x\right\|_{E}\leqslant 1\right\\}\;.$

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

$\displaystyle\left\langle\ell(x),y\right\rangle_{E}=\left\langle
x,\ell^{T}(y)\right\rangle_{F}\qquad\text{for all $(x,y)\in F\times E$}$ (57)

and, for convenience, we also denote

$\displaystyle\ell^{\nabla}:=\ell^{T}\circ\nabla A\circ\ell\;.$ (58)

Here the function $A$ is the one appearing in 1 an exponential family
distribution. Recall that $A$ is infinitely differentiable on
$\mathrm{Int}E_{0}$ so that by (C1), $\ell^{\nabla}$ is well-defined on
$\Upsilon(\tilde{{\mathsf{T}}})$ for $\tau$-a.e. $\ell$ with
$\tau\in\mathrm{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 $\mathsf{Y}\times\mathsf{Z}$,

$\displaystyle\mathbb{E}_{\vartheta,\tau}\left[g(Y,L)\right]=\int
g(y,\ell)\,k^{(o)}(\ell\circ\Upsilon(\vartheta),y)\;\nu(\mathrm{d}y)\tau(\mathrm{d}\ell)\;.$
(59)

We now introduce, for any
$(\vartheta_{0},\tau_{0})\in\tilde{{\mathsf{T}}}\times\mathrm{M}$, two
auxiliary functions $\tilde{\mathrm{m}}_{\vartheta_{0},\tau_{0}}$ and
$\tilde{\mathrm{s}}_{\vartheta_{0}\tau_{0}}$, both defined from $\mathsf{Y}$
to $\mathbb{R}_{+}$, such that the two following conditions hold.

  1. (C2)

For all $\vartheta_{0}\in\tilde{{\mathsf{T}}}$ and $\tau_{0}\in\mathrm{M}$, we
have

$\mathbb{E}_{\vartheta_{0},\tau_{0}}\left[\left.\left\|L\right\|_{\mathrm{op}}\,\right|\,Y\right]\leqslant\tilde{\mathrm{m}}_{\vartheta_{0},\tau_{0}}(Y)\quad\mathbb{P}_{\vartheta_{0},\tau_{0}}-\mathrm{a.s.}$

  2. (C3)

For all $\vartheta_{0}\in\tilde{{\mathsf{T}}}$ and $\tau_{0}\in\mathrm{M}$, we
have

$\forall\xi\in\Upsilon(\tilde{{\mathsf{T}}})\,,\,\exists\epsilon,C>0\,,\,\mathbb{E}_{\vartheta_{0},\tau_{0}}\left[\left.\mathrm{e}^{\epsilon\left\|L^{\nabla}(\xi)\right\|_{F}}\,\right|\,Y\right]\leqslant
C\,\mathrm{e}^{\tilde{\mathrm{s}}_{\vartheta_{0},\tau_{0}}(Y)}\quad\mathbb{P}_{\vartheta_{0},\tau_{0}}-\mathrm{a.s.}$

Assumptions (C2) and (C3) above are made to ensure that the integrals used to
solve the argmax problem from 9 will be well-defined.

We finally present our last assumption, which corresponds to some sort of
identifiability assumption (see 10 in Section C.4 for details). It is used to
obtain the uniqueness for the argmax problem in 5.

  1. (C4)

For all $\xi\neq\xi^{\prime}\in\Upsilon(\tilde{{\mathsf{T}}})$ and
$\tau_{0}\in\mathrm{M}$, we have
$\tau_{0}\left(\left\\{\ell\in\mathsf{Z}\;:\;\ell(\xi)\neq\ell(\xi^{\prime})\right\\}\right)>0$.

We can now state the following result, whose proof can be found in Section
C.4.

######  Theorem 5.

Consider an ELM family as in 2 and let
$\theta_{0}=(\vartheta_{0},\tau_{0})\in\tilde{{\mathsf{T}}}\times\mathrm{M}$.
Assume (C1)–(C4) hold. Let $\check{\varphi}$ be a probability density function
w.r.t. $\nu$ such that

$\int_{\mathsf{Y}}\check{\varphi}(y)\,\left(\left\|S(y)\right\|_{E}\,\tilde{\mathrm{m}}_{\vartheta_{0},\tau_{0}}(y)+\mathrm{e}^{\tilde{\mathrm{s}}_{\vartheta_{0},\tau_{0}}(y)}\right)\;\nu(\mathrm{d}y)<\infty\;,$
(60)

and define $\varphi:\mathsf{Z}\times\mathsf{Y}\to\mathbb{R}_{+}$ by (54). Let
us further set for all
$(\ell,y,\xi)\in\mathsf{Z}\times\mathsf{Y}\times\Upsilon(\tilde{{\mathsf{T}}})$
such that $\ell\circ\Upsilon(\vartheta_{0})\in\mathrm{Int}E_{0}$,

$\displaystyle\tilde{\varphi}(\ell)=\int_{\mathsf{Y}}\varphi(\ell,y)\;\nu(\mathrm{d}y)\;,$
(61)
$\displaystyle\mathrm{w}^{\nabla}_{\tau_{0}}(\xi)=\int_{\mathsf{Z}}\tilde{\varphi}(\ell)\,\ell^{\nabla}(\xi)\;\tau_{0}(\mathrm{d}\ell)\;.$
(62)

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

$\displaystyle\mathrm{w}^{\nabla}_{\tau_{0}}\circ\Upsilon(\vartheta)=\mathbb{E}_{\vartheta,\tau_{0}}\left[L^{T}\circ
S(Y)\,\tilde{\varphi}(L)\right]\;.$ (63)

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

$\mathop{\mathrm{argmax}}_{\xi\in\Upsilon(\tilde{{\mathsf{T}}})}\int_{\mathsf{Y}}\left(\int_{\mathsf{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}(\mathrm{d}\ell)\right)\nu(\mathrm{d}y)\;,$
(64)

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(\mathrm{d}y)\tau_{0}(\mathrm{d}\ell)+\frac{b}{1+b}\,\mathrm{w}^{\nabla}_{\tau_{0}}(\xi_{0})$
(65)

belongs to
$\mathrm{w}^{\nabla}_{\tau_{0}}\circ\Upsilon(\tilde{{\mathsf{T}}})$, in which
case (64) has a unique solution $\xi^{*}$ defined by

$\mathrm{w}^{\nabla}_{\tau_{0}}(\xi^{*})=\mathbf{w}^{*}\;.$ (66)

Observe that by setting $\xi^{*}=\Upsilon(\vartheta^{*})$ and
$\xi_{0}=\Upsilon(\vartheta_{0})$ as well as using (63), the equation (66) can
be interpreted as

$\mathbb{E}_{\vartheta^{*},\tau_{0}}\left[L^{T}\circ
S(Y)\,\tilde{\varphi}(L)\right]=\mathbb{E}_{\vartheta_{0},\tilde{\tau}_{0}}\left[L^{T}\circ
S(\tilde{Y})\right]\;,$ (67)

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 $\mathsf{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\left[\left.y\,\right|\,\ell\right]=\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 (67) as
well, the latter reads

$\mathbb{E}_{\vartheta^{*},\tilde{\tau}_{0}}\left[L^{T}\circ
S(Y)\right]=\mathbb{E}_{\vartheta_{0},\tilde{\tau}_{0}}\left[L^{T}\circ
S(\tilde{Y})\right]\;.$ (68)

Hence, 5 along with (68) can be used to solve the argmax problem (56) 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 (60) on
$\check{\varphi}$ holds. This can notably be done when $\check{\varphi}$ takes
one of the two forms listed in 1 (see 8 in Section 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 $\mathsf{Y}=\mathbb{R}^{d}$, with
$d\geqslant 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)\;\mathrm{d}y<\infty\;.$
(69)

We let $J\geqslant 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{\mathsf{T}}=\mathbb{R}^{d}\times\mathcal{M}_{>0}(d)\times(0,\infty)$,
defined by

$k(\theta,y)=\frac{\Gamma\left((a+d)/2\right)}{\Gamma\left(a/2\right)(a\pi)^{d/2}\left|\Sigma\right|^{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\geqslant 1}$ valued in $\mathcal{Q}$ where
$\mu_{n}=\mu_{{\boldsymbol{\lambda}_{n}},\Theta_{n}}$ with
${\boldsymbol{\lambda}_{n}}\in\mathcal{S}_{J}^{+}$ and
$\Theta_{n}=(m_{j,n},\Sigma_{j,n},a_{j,n})_{j=1\dots J}\in{\mathsf{T}}^{J}$
such that, for all $n\geqslant 1$,

  1. (i)

given $(\eta_{n})_{n\geqslant 1}$ valued in $(0,1]$, $(\kappa_{n})_{n\geqslant
1}$ valued in $(-\infty,0]$ and
${\boldsymbol{\lambda}_{1}}\in\mathcal{S}_{J}^{+}$,
$({\boldsymbol{\lambda}_{n}})_{n\geqslant 2}$ is defined by the iterative
formula (23);

  2. (ii)

given $(\gamma_{j,n})_{j=1,\dots,J,n\geqslant 1}$ valued in $(0,1]$, and
$\Theta_{1}\in{\mathsf{T}}^{J}$, $(\Theta_{n})_{n\geqslant 2}$ is defined by
the iterative formula:

$\displaystyle
m_{j,n+1}=\left(\int_{\mathbb{R}^{d}\times\mathbb{R}_{+}}z\,\psi_{j,n}(z,y)\;\mathrm{d}y\,\check{\tau}_{a_{j,n}}(\mathrm{d}z)\right)^{-1}\,\int_{\mathbb{R}^{d}\times\mathbb{R}_{+}}y\,z\,\psi_{j,n}(z,y)\;\mathrm{d}y\,\check{\tau}_{a_{j,n}}(\mathrm{d}z)\;,$
(70)
$\displaystyle\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)\;\mathrm{d}y\,\check{\tau}_{a_{j,n}}(\mathrm{d}z)\;,$
(71) $\displaystyle
a_{j,n+1}=2\,\kappa^{-1}\left(\int_{\mathbb{R}^{d}\times\mathbb{R}_{+}}(z-\ln(z))\,\varphi_{j,n}(z,y)\;\mathrm{d}y\,\check{\tau}_{a_{j,n}}(\mathrm{d}z)\right)\;,$
(72)

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^{\prime}(x)/\Gamma(x)$, we set

$\displaystyle\varphi_{j,n}(z,y)=\check{\varphi}_{j,n}^{(\alpha)}(y)\frac{\mathrm{e}^{\frac{-z}{2}(y-m_{j,n})^{T}\Sigma_{j,n}^{-1}(y-m_{j,n})}}{(2\pi/z)^{d/2}\left|\Sigma_{j,n}\right|^{1/2}k(\theta_{j,n},y)}\;,$
(73)
$\displaystyle\psi_{j,n}(z,y)=\gamma_{j,n}\,\varphi_{j,n}(z,y)+(1-\gamma_{j,n})\,\,\frac{\mathrm{e}^{\frac{-z}{2}(y-m_{j,n})^{T}\Sigma_{j,n}^{-1}(y-m_{j,n})}}{(2\pi/z)^{d/2}\left|\Sigma_{j,n}\right|^{1/2}}\,\left(\int\varphi_{j,n}(z,y)\;\mathrm{d}y\right)\;,$
(74)
$\displaystyle\check{\tau}_{a}(\mathrm{d}z)=\mathbbm{1}_{\mathbb{R}_{+}}(z)\,\frac{(a/2)^{a/2}}{\Gamma(a/2)}z^{a/2-1}\mathrm{e}^{-za}\;\mathrm{d}z\quad\text{for
any $a>0$}\;,$ (75)

Observe that $\check{\tau}_{a}$ in (75) is the $\chi^{2}$ distribution with
$a$ degrees of freedom and that $\kappa^{-1}$ in (72) is well-defined as a
consequence of 13.

We then have the following result, whose proof is postponed to Section C.4.

######  Corollary 10.

In the setting of 5, if $\Psi_{\alpha}(\mu_{1}k)<\infty$, then, at any time
$n\geqslant 1$, we have
$\Psi_{\alpha}(\mu_{n+1}k)\leqslant\Psi_{\alpha}(\mu_{n}k)$.

Notice that all the update formulas established so far involve intractable
integrals. In the next section, we focus on the particular case of Gaussian
Mixture Models seen in Examples 3 and 4 and we investigate how our framework
can be used to derive tractable algorithms for Gaussian Mixture Model
optimisation.

## 7 Gaussian Mixture Models (GMMs) optimisation

We consider the case of $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{\mathsf{T}}$ denotes the mean and covariance
matrix of the $j$-th Gaussian component density. Recall that, by 2, a valid
update formula at time $n\geqslant 1$ for the mixture weights
$({\boldsymbol{\lambda}_{n}},\Theta_{n})_{n\geqslant 1}$ is

$\displaystyle\lambda_{j,n+1}$
$\displaystyle=\frac{\lambda_{j,n}\left[\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}{\sum_{\ell=1}^{J}\lambda_{\ell,n}\left[\int_{\mathsf{Y}}\varphi_{\ell,n}^{(\alpha)}(y)\nu(\mathrm{d}y)+(\alpha-1)\kappa_{n}\right]^{\eta_{n}}}\;,\quad
j=1\ldots J\;.$ (76)

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

$\displaystyle\forall j=1\ldots J,\quad m_{j,n+1}$
$\displaystyle=(1-\gamma_{j,n})m_{j,n}+\gamma_{j,n}\widehat{m}_{j,n}$ (77)
$\displaystyle\Sigma_{j,n+1}$
$\displaystyle=(1-\gamma_{j,n})\Sigma_{j,n}+\gamma_{j,n}\widehat{\Sigma}_{j,n}+\gamma_{j,n}(1-\gamma_{j,n})\left(\widehat{m}_{j,n}-m_{j,n}\right)\left(\widehat{m}_{j,n}-m_{j,n}\right)^{T}$

where

$\displaystyle\widehat{m}_{j,n}=\int_{\mathsf{Y}}y\,\check{\varphi}_{j,n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)\quad\mbox{and}\quad\widehat{\Sigma}_{j,n}=\int_{\mathsf{Y}}yy^{T}\,\check{\varphi}_{j,n}^{(\alpha)}(y)\;\nu(\mathrm{d}y)-\widehat{m}_{j,n}\widehat{m}_{j,n}^{T}$

with
$\check{\varphi}_{j,n}^{(\alpha)}=\varphi_{j,n}^{(\alpha)}/\int\varphi_{j,n}^{(\alpha)}\nu(\mathrm{d}y)$
and $\gamma_{j,n}\in(0,1]$. In addition, by virtue of 4, another possible
update for the means corresponds to setting:

$\displaystyle m_{j,n+1}$
$\displaystyle=m_{j,n}+\gamma_{j,n}\frac{\int_{\mathsf{Y}}\lambda_{j,n}\varphi_{j,n}^{(\alpha)}(y)(y-m_{j,n})\nu(\mathrm{d}y)}{\int_{\mathsf{Y}}{\mu_{n}k(y)^{\alpha}p(y)^{1-\alpha}}\nu(\mathrm{d}y)}\;,\quad
j=1\ldots J\;$ (78)

in lieu of the means update written in (77) and with $\gamma_{j,n}$ still in
$(0,1]$. Note that since the update (78) is linked to Gradient Descent steps
for Rényi’s $\alpha$-divergence minimisation (see 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 (77) will be referred to as the
Maximisation Gradient (MG) update.

At this stage, one may observe that the updates (76), (77) and (78) all
involve intractable integrals. We can bypass this difficulty by resorting to
approximate update rules and by taking advantage of the fact that the
integrals appearing in those update formulas are of the form

$\displaystyle\int_{\mathsf{Y}}\varphi_{j,n}^{(\alpha)}(y)g(y)\nu(\mathrm{d}y)\;,$
(79)

where $g$ is a well-chosen function defined on $\mathsf{Y}$ and where the
definition of the function $\varphi_{j,n}^{(\alpha)}$ in (16) is recalled
below: for all $y\in\mathsf{Y}$,

$\displaystyle\varphi_{j,n}^{(\alpha)}(y)=k(\theta_{j,n},y)\left(\frac{\mu_{n}k(y)}{p(y)}\right)^{\alpha-1}\;.$

From there, many choices are possible and for simplicity we restrict ourselves
to performing typical Importance Sampling estimation. Looking at the
definition of $\varphi_{j,n}^{(\alpha)}$, a first idea could be to sample $M$
i.i.d random variables $(Y_{m,n}^{(j)})_{1\leqslant m\leqslant M}$ according
to $k(\theta_{j,n},\cdot)$ for $j=1\ldots J$ and to use the unbiased estimator
of (79) given by

$\displaystyle\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)})\;.$

However, this sampler depends on $j$ and this option thus 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\geqslant n}$ that do not depend on $j$ at
time $n$. More specifically, we discuss two possibilities.

  1. (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\leqslant m\leqslant M}$
according to the best approximation of $p$ we have a time $n$, that is
$q_{n}=\mu_{n}k$. In that case, the samples will be 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 (79) is given by:

$\displaystyle\frac{1}{M}\sum_{m=1}^{M}\hat{\varphi}_{j,n}^{(\alpha)}(Y_{m,n})g(Y_{m,n})\;,$
(80)

where
$\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}$
for all $y\in\mathsf{Y}$ and all $j=1\ldots J$.

  2. (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 (80). 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
${\boldsymbol{\lambda}_{n}}$.

######  Remark 3.

On a side note, approximated update rules for the GMM case can also be
obtained by sampling according to the sampler
$q_{n}=\mathcal{N}(\cdot;\boldsymbol{0}_{d},\boldsymbol{I}_{d})$, where
$\boldsymbol{0}_{d}$ is the null vector of dimension $d$. Indeed, those can be
deduced by applying the reparameterisation trick: for all $m=1\ldots 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
(79) is then given by:

$\displaystyle\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)})\;.$

Contrary to (i) and (ii), this choice of sampler relies on the existence of a
reparameterisation. 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.

Note that our framework 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 below.

At iteration $n$,

  1. 1.

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

  2. 2.

For all $j=1\ldots J$,

    1. (a)

Compute

$\displaystyle\widehat{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})}\leavevmode\nobreak\
\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\ \leavevmode\nobreak\
\widehat{\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})}-\widehat{m}_{j,n}\widehat{m}_{j,n}^{T}$

    2. (b)

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

$\displaystyle(MG)\quad m_{j,n+1}$
$\displaystyle=\left(1-\gamma_{n}\right)m_{j,n}+\gamma_{n}\widehat{m}_{j,n}$
$\displaystyle(RGD)\quad m_{j,n+1}$
$\displaystyle=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

$\displaystyle\lambda_{j,n+1}$
$\displaystyle=\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}}}$
$\displaystyle\Sigma_{j,n+1}$
$\displaystyle=(1-\gamma_{n})\Sigma_{j,n}+\gamma_{n}\widehat{\Sigma}_{j,n}+\gamma_{n}(1-\gamma_{n})\left(\widehat{m}_{j,n}-m_{j,n}\right)\left(\widehat{m}_{j,n}-m_{j,n}\right)^{T}$

Algorithm 3 GMMs optimisation with the IS-n/IS-unif sampler

######  Remark 4.

Computing the estimator (80) 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,\mathscr{D})$ often consider large-scale datasets
$\mathscr{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,\mathscr{D})$, as detailed in [19,
Section 4.3].

Now that we have presented an example of practical implementations of our
algorithms, we move on to presenting our numerical experiments.

## 8 Numerical experiments

In our numerical experiments, we are interested in seeing how the choice of
our hyperparameters influences the convergence of Algorithm 3. We will
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.

Letting $\boldsymbol{u_{d}}$ be the $d$-dimensional vector whose coordinates
are all equal to $1$, $\boldsymbol{I_{d}}$ be the identity matrix and $c$ be a
positive constant (we set $c=2$), we will consider three multimodal examples:

  1. (i)

Equally-weighted Gaussian Mixture Model. The target $p$ is a mixture density
of two equally-weighted $d$-dimensional Gaussian distributions multiplied by
the positive constant $c$ such that

$p(y)=c\times\left[0.5\mathcal{N}({y};-2\boldsymbol{u_{d}},\boldsymbol{I_{d}})+0.5\mathcal{N}({y};2\boldsymbol{u_{d}},\boldsymbol{I_{d}})\right]\;.$

  2. (ii)

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

$p(y)=c\times\left[0.35\mathcal{N}({y};-2\boldsymbol{u_{d}},\boldsymbol{I_{d}})+0.25\mathcal{N}({y};2\boldsymbol{u_{d}},\boldsymbol{I_{d}})+0.4\mathcal{N}({y};\boldsymbol{u_{d}},\boldsymbol{I_{d}})\right]$

  3. (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 the positive constant
$c$ such that