A Dynamical Systems Approach for Convergence of the Bayesian EM Algorithm

1 Introduction

This work is inspired by the recent stream of papers importing ideas from (dynamical) systems and control (S&C) theory into optimization [Elia2011, Su2014, Lessard2016, Wibisono2016, Fazlyab2017b, Scieur2017, Franca2018, Taylor2018, Wilson2018, Franca2019, Orvieto2019]. While a significant volume of these optimization-based papers have been published at machine learning (ML) venues, only a few have been explicitly dedicated to addressing concrete ML problems, applications, or algorithms [Plumbley1995, Pequito2011UnsupervisedLO, Zhu2018AnOC, Aquilanti2019, Liu2019]. Furthermore, only a small subset of this emerging topic of research has focused directly on discrete-time analysis that is the direct result of discretizations of an underlying continuous-time version of the algorithms [Lessard2016, Fazlyab2018, Fazlyab2018b, Taylor2018, romero2019].

The aforementioned effort to bring ideas from S&C into optimization is largely done via Lyapunov-based tools, which are most popular for nonlinear systems [Khalil2001]. Simply put, Lyapunov functions may be seen as abstract surrogates of energy in a dynamical system. If the function is monotonically decreasing over time, then some form of stability must be present. Likewise, if persistently strictly increasing, then instability is inevitable. There exist several close relationships and parallels between notions of stability of dynamical systems and convergence of optimization algorithms. However, to the best of the authors’ knowledge, the current literature lacks a comprehensive summary of these relationships, with the closest work that we are aware of being are [Schropp1995, Lessard2016, Fazlyab2018, Taylor2018, Wilson2018].

Furthermore, most literature in S&C regarding Lyapunov stability theory deals exclusively with continuous-time systems, and (dynamical) stability properties are typically not preserved after discretization [Stuart1998]. However, a vast amount of general Lyapunov-based results possess a discrete-time equivalent [Kalman1959, Byrnes1994, Agarwal2000, Jiang2001, Jiang2002, Lakshmikantham2002, Kellett2004, Kellett2005, Geiselhart2016, Bof2018]. These discrete-time equivalent Lyapunov-based stability results have largely not been adequately leveraged for convergence analysis of optimization algorithms with ML applications.

Further challenges arise when it comes to putting Lyapunov’s ideas to practice, particularly in the context of O&ML. Most notably, finding a suitable Lyapunov function, often out of a set of candidates, may prove extraordinarily difficult. For instance, for physical engineering systems, Lyapunov functions are often constructed from first principles as some meaningful measure of energy. In optimization, the situation is somewhat similar in the sense that a suitable Lyapunov function may often be constructed as a surrogate of the objective function. There exist some general Lyapunov function construction methods, but they steep come with natural limits regarding their scope [Li2014]. Outside of these scenarios, Lyapunov functions are traditionally constructed rather empirically, on a case-by-case basis. Furthermore, while there exist in the literature certain reciprocal existence results of Lyapunov functions for stable systems, as well as Lyapunov-based instability results, it should be noted that not being able to find or construct a Lyapunov function is not enough ground to assert instability.

Contributions

In this paper, we conduct a S&C-based analysis of the convergence of a widely popular algorithm used for incomplete-data estimation and unsupervised learning – the expectation-maximization (EM) algorithm. More precisely, we will focus on the Bayesian variant of the EM algorithm originally proposed in

[Dempster1977], which we refer to as the MAP-EM algorithm, since it is used for maximum a posteriori (MAP) estimation [Figueiredo2004]. We will leverage notions from discrete-time Lyapunov stability theory to study the convergence of MAP-EM, and in the process provide exclusive insights on the robustness of our derived conditions for stability (asymptotic or otherwise), and thus convergence guarantees. Furthermore, it should be noted that these results could likely not have been constructed without our approach. The primary contributions of this paper are:

[noitemsep,topsep=0pt]
Using Lyapunov stability theory, we show that, under certain regularity assumptions, the MAP-EM algorithm is locally convergent to strict local maxima of the incomplete-data posterior, including the MAP estimator.
We provide conditions to establish Lyapunov exponential stability, which translates to a linear convergence rate of EM, or even a quadratic one under stronger conditions.
We prove that under certain conditions (most notably, unimodality of the incomplete-data posterior), the EM algorithm is globally convergent.

With this, we argue for the possibility of extending our S&C-based framework to discover robust stability conditions and novel convergence guarantees of alternative iterative optimization algorithms used in machine learning.

2 Background: MAP-EM Algorithm

Let $θ$ be an unknown parameter of interest that we seek to infer from an idealized unobservable dataset $x$ , via the statistical model $p (x | θ)$ and prior $p (θ)$ . Given that $x$

is not directly observable, suppose another random variable

y

, which may be seen as an incomplete version of

x

, is observable. In this definition,

x = (y, z)

, with

z

seen as missing data. For this reason,

x

is typically referred to as the complete dataset. More generally, we could have

x = g (y, z)

for some

g (\cdot)

, with

y

observable and

z

hidden, or, simplistically

y = h (x)

for some

h (\cdot)

. In practice, there could even be no explicit missing data or no relationship between

x

and

y

in terms of transformations. Instead the only relationship could be the Markov condition

θ \to x \to y

[Gupta2011], meaning that

y

is conditionally independent of

θ

subject to

x

, i.e.

p (y | x, θ) = p (y | x)

We adopt the notation that $x, y, θ$

are all (absolutely) continuous random variables, but, in reality, only

θ | y

is required to be so, with

x

and

y

being allowed to be discrete, continuous, or of mixed type. For ease of notation, we use

d x, d y, d θ

to refer, respectively, to integration with respect to (w.r.t.) implicit

σ

-finite measures that dominate the probability distributions of

x, y, θ

, or appropriate conditionals of these, and which coincide with the construction of the respective densities via Radon-Nikodym derivatives. We aim to estimate

θ

from

y

via the maximum a posteriori (MAP) estimator:

{^θ}_{MAP} ≜ {a r g m a x}_{θ} log p (θ | y)

, where the maximization is taken over the entire parameter space. In some situations a global maximizer may not exist, and we need to be content with (or even give preference to) a “good” local maximizer [Figueiredo2004]. The mapping

θ \mapsto log p (θ | y)

is typically referred to as the incomplete-data log-likelihood function, whereas

θ \mapsto log p (θ | x)

the complete-data log-likelihood function.

Under the incomplete-data framework, a popular approach to compute the MAP estimator is through the expectation-maximization (EM) algorithm (also known as MAP-EM in our Bayesian framework), whose iterations are

{^θ}_{k + 1} = a r g m a x θ Q (θ, {^θ}_{k}), where, Q (θ,^θ) ≜ E_{x \sim p (x | y,^θ)} [log p (θ | x)]

(1)

from a given initial estimate ${^θ}_{0}$ , where $Q (θ,^θ)$ denotes the expected complete-data log-posterior, conditional to the observed data $y$ and current parameter estimate $^θ$ . The maximum in (1) is taken within the entire parameter space. Eventually, as we will demonstrate later, EM is (in general) a local (greedy) search method w.r.t. the actual objective function – the (incomplete-data) log-posterior ¹¹1The underlying assumption of the EM algorithm is that (1) is easy to globally maximize in $θ$ . If it can’t be exactly and globally maximized, but instead if we can find some ${^θ}_{k + 1}$ such that $Q ({^θ}_{k + 1}, {^θ}_{k}) > Q ({^θ}_{k}, {^θ}_{k}), \forall θ,$ (potentially excluding the case when ${^θ}_{k}$ is already fixed point of EM), then any variant of EM that settles for such a sequence ${^θk}k∈Z+$ of iterates, where $Z_{+} = {0, 1, 2, \dots}$ , is known as a generalized EM (GEM) algorithm..

2.1 An Information-Theoretic Perspective

As we dive into the information theoretic perspective, recall that $D_{K L} (p ∥ q) ≜ \int_{X} p (x) log (\frac{}{p (x)} q (x)) d x$

denotes the Kullback-Leibler (KL) divergence between the probability density functions

p (x)

and

q (x)

w.r.t. the same

σ

-finite dominating measure over

X

implicitly denoted via

d x

. The (random) variables of interest satisfies the following assumption

Assumption 1.

The random variables $x, y, θ$ satisfy the Markovian condition $θ \to x \to y$ .

Now, the Markov Assumption 1, i.e. $p (y | x, θ) = p (y | x)$ , allows us to restate the $Q$ -function and subsequently the EM algorithm in information-theoretic terms, leading to the following proposition.

Proposition 1.

Under Assumption 1, we have

\begin{matrix} Q (θ,^θ) = log p (θ | y) - D_{K L} [p_{X} (\cdot | y,^θ) ∥ p_{X} (\cdot | y, θ)] + terms that do not depend on θ, \end{matrix}

(2)

and thus, the MAP-EM algorithm (1) can be rewritten as

^θk+1=argmaxθ{logp(θ|y)−d(θ,^θk)},

(3)

where $d (θ,^θ) ≜ D_{K L} [p_{X} (\cdot | y,^θ) ∥ p_{X} (\cdot | y, θ)]$ and $p_{X} (x | y, θ) ≜ p (x | y, θ)$ .

Since the sole purpose of the $Q$ -function is to be maximized at the M-step, we thus redefine it as $Q (θ,^θ) ≜ log p (θ | y) - d (θ,^θ)$ by dropping the terms that do not depend on $θ$ . Notice that the EM algorithm can be recognized as a (generalized) proximal point algorithm (PPA)

^θk=argminθ{ℓ(θ)+βkd(θ,^θ)},

(4)

with a fixed step size $β_{k} = 1$ , where the function we seek to minimize $ℓ (θ) ≜ - log p (θ | y)$ . Furthermore, $d (\cdot, \cdot)$ may be seen as a regularizing premetric²²2It is well known that $d (θ,^θ) \geq 0$ , with equality corresponding almost exclusively to $θ =^θ$ . More precisely, consider two probability distributions $P$ and $Q$ over the same measurable space $(S, F)$ , both dominated over some $σ$ -finite measure $ν$ in the space $X$ . With a slight overload of notation, let $D_{K L} (P ∥ Q) ≜ E_{P} [log (\frac{d P}{d Q})] = \int_{S} log (\frac{d P}{d Q}) d P = \int_{X} p log (\frac{p}{q}) d ν$ , where $p = \frac{d P}{d ν}$ and $q = \frac{d Q}{d ν}$ . Then, $D_{K L} (P ∥ Q) \geq 0$ with equality if and only if $\frac{d P}{d Q} = 1$ ( $P$ -a.s.), which is also equivalent to $p = q$ ( $P$ -a.s.). Despite this, the KL divergence is non-symmetric and does not satisfy the triangle inequality, in general, and thus it is not a metric.. This perspective was explored in detail in [Chretien2000, Figueiredo2004], which served as inspiration for this work.

Notice also that, if the prior $p (θ)$

is “flat”, meaning that we assign the (typically) non-informative (and potentially degenerate) prior given by a uniform distribution over the parameter space, then

ℓ (θ) \propto - log p (y | θ)

, with

θ \mapsto log p (y | θ)

naturally denoting the incomplete-data log-likelihood function. Therefore, the Bayesian EM algorithm (MAP-EM) generalizes the traditional EM algorithm, and for this reason, from this point on we will largely refer to MAP-EM as simply EM.

3 The Dynamical Systems Approach

We now interpret the EM algorithm as a (discrete-time and time-invariant) nonlinear state-space dynamical system,

{^θ}_{k + 1}

= F ({^θ}_{k}), \forall k \in Z_{+}, % w i t h, F (^θ) ≜ a r g m a x θ Q (θ,^θ) .

(5)

In the language of dynamical systems, ${^θ}_{k}$ is known as the state of the system (5) at time step $k$ (usually denoted as $x_{k}$

, much like in reinforcement learning, but for the sake of notational consistency with EM we opted for

{^θ}_{k})

, and, in particular,

{^θ}_{0}

is called the initial state. The space of points from which the state can take values is known as the state space. In our case, the parameter space and the state space coincide. The sequence

{^θ0,^θ1,…}

is often called a trajectory starting from initial state

{^θ}_{0}

. The function

F

represents the dynamics of the system. The following assumption ensures that

F (^θ)

is uniquely defined.

Assumption 2.

$Q (\cdot,^θ)$ has a unique global maximizer for each $^θ$ .

In other words, the complete-data log-posterior is expected to have a unique global maximizer, which in principle does not prevent the incomplete-data log-posterior from having multiple global maxima or from being unbounded, such as in the cases of GMMs with unknown covariance matrices.

3.1 Equilibrium in Dynamical Systems

We say that a point ${^θ}^{⋆}$ in the state space is an equilibrium of the dynamical system (5) if ${^θ}_{0} = {^θ}^{⋆}$ implies that ${^θ}_{k} = {^θ}^{⋆}$ for every $k \in Z_{+}$ . In other words, ${^θ}^{⋆}$ is an equilibrium point if and only if ${^θ}^{⋆}$ is a fixed point of $F$ . This implies that, $F ({^θ}^{⋆}) = {^θ}^{⋆}$ , the equilibria of the dynamical system representation of EM naturally coincide with its fixed points.

Recall that limit points of the EM algorithm (1) consist of points ${^θ}^{⋆}$ for which there exists some ${^θ}_{0}$ such that ${^θ}_{k} \to {^θ}^{⋆}$ as $k \to \infty$ . However, EM need not be locally convergent near limit points of its dynamical system representation, which requires that ${^θ}_{k} \to {^θ}^{⋆}$ as $k \to \infty$ for every ${^θ}_{0}$ in a small enough neighborhood of ${^θ}^{⋆}$ . Such points are known as (locally) attractive in dynamical systems theory. In order to establish a key relationship between limit points of EM and equilibria of its dynamical systems representation, $F$ needs to be continuous as stated in Assumption 3.

Assumption 3.

$Q (\cdot)$ is continuous.

This follows, for instance, if $θ \mapsto p (θ | y)$ is continuous and $x$ conditional to $y$ and $θ$ has a finite support, which is the case for GMM clustering.

Lemma 1.

Under Assumptions 1–3, $F$ is continuous.

Referring the readers to supplementary section for Lemma 1’s proof, we move ahead to establish a key relationship between limit points of EM and its dynamical systems representation.

Proposition 2.

Under Assumptions 1–3, any limit point of the EM algorithm (1) is also a fixed point of $F$ given by (5), and thus an equilibrium of EM’s dynamical system representation (5).

The reciprocal is not true, however, which is known to occur for unstable equilibria of nonlinear systems. Therefore, we now focus on another key concept in dynamical systems theory – that of (Lyapunov) stability – and proceed to study its relationship with convergence of the EM algorithm.

3.2 Lyapunov Stability

We say that ${^θ}^{⋆}$ in the parameter space is a stable point of the system (5) if the trajectory ${^θ0,^θ1,^θ2,…}$ is arbitrarily close to ${^θ}^{⋆}$ , provided that $∥ {^θ}_{0} - {^θ}^{⋆} ∥ > 0$ is sufficiently small. In means that, if and only, for any $ε > 0$ , there exists some $δ > 0$ such that for every ${^θ}_{0}$ satisfying $∥ {^θ}_{0} - {^θ}^{⋆} ∥ \leq δ$ , we have $∥ {^θ}_{k} - {^θ}^{⋆} ∥ \leq ε$ for every $k \in Z_{+}$ . If, in addition, there exists some $δ > 0$ small enough such that, for every ${^θ}_{0}$ satisfying $∥ {^θ}_{0} - {^θ}^{⋆} ∥ \leq δ$ , we have ${^θ}_{k} \to {^θ}^{⋆}$ as $k \to \infty$ , then we say ${^θ}^{⋆}$ is a (locally) asymptotically stable point of the system. In other words, asymptotically stable points are simply stable and attractive points of the system. If the attractiveness is global, meaning that $δ > 0$ can be made arbitrarily large, then we say that ${^θ}^{⋆}$ is globally asymptotically stable. Finally, if in addition there exist $δ, c > 0$ and $ρ \in (0, 1)$ such that, for every ${^θ}_{0}$ satisfying $∥ {^θ}_{0} - {^θ}^{⋆} ∥ \leq δ$ , we have $∥ {^θ}_{k} - {^θ}^{⋆} ∥ \leq c ρ^{k}$ for every $k \in Z_{+}$ , then we say that ${^θ}^{⋆}$ is a (locally) exponentially stable point of the system. Global exponential stability holds if $δ > 0$ can be made arbitrarily large.

Lemma 2.

Consider the dynamical system (5) with an arbitrary continuous function $F$ . Then, every stable point is also an equilibrium.

In particular, we can now see that, under Assumptions 1–3, every stable point of the dynamical system that represents EM is also a fixed point of EM. Furthermore, it should be clear that local maxima of the incomplete-data log-posterior (that happen to be asymptotically stable) must be locally convergent points of EM. On the other hand, exponentially stable local maxima lead EM to attain a $ρ$ -linear convergence rate, where $ρ$ originates from the definition of exponential stability. In particular, if the cost function $ℓ (^θ)$ is Lipschitz continuous, then clearly $ℓ ({^θ}_{k}) - ℓ ({^θ}^{⋆}) = O (ρ^{k})$ . The the relationships between all the aforementioned concepts is summarized in the supplementary section B.

To establish the different notions of stability, we use the ideas proposed by Lyapunov, which we summarize in Lemma 3 in what we refer to as the “Lyapunov theorem” (actually a collective of results). Before stating the Lyapunov theorem, let us introduce some convenient terminology borrowed from nonlinear systems theory. We say that a function $V (^θ)$ is:

[noitemsep,topsep=-5pt]
(locally) positive semidefinite w.r.t. ${^θ}^{⋆}$ if $V ({^θ}^{⋆}) = 0$ and $V (^θ) \geq 0$ for $^θ$ near ${^θ}^{⋆}$ ;
(locally) positive definite w.r.t ${^θ}^{⋆}$ if $V (^θ) \geq 0$ with equality if and only if $^θ =^θ$ , for every $^θ$ near ${^θ}^{⋆}$ ;
(locally) negative semidefinite (respectively, negative definite) if $- V$ is positive semidefinite (respectively, positive definite);
radially unbounded if the state space is an entire Euclidean space $R^{p}$ and $V (^θ) \to + \infty$ as $∥^θ ∥ \to \infty$ .

Global variants of definitions 1–3 hold if $^θ$ can be picked anywhere in the state space. Finally, we say that a scalar function $α : [0, \infty) \to [0, \infty)$ is of class $K$ if $α (0) = 0$ and if it is continuous and strictly increasing.

We are now ready to state Lyapunov’s theorem for generic discrete-time systems of the form (5). These results, and Lyapunov-based results alike, are of crucial importance in practice for showing the different forms of stability of nonlinear systems.

Lemma 3 (Lyapunov theorem).

Consider the dynamical system (5) with an arbitrary continuous function $F$ , and let ${^θ}^{⋆}$ be an arbitrary point in the state space. Let $V (^θ)$ be a continuous function and $Δ V ≜ V \circ F - V$ . Consider the following assertions about these functions:

[noitemsep,topsep=0pt]
$V$ is positive definite w.r.t. ${^θ}^{⋆}$ ;
$Δ V$ is negative semidefinite w.r.t. ${^θ}^{⋆}$ ;
$Δ V$ is negative definite w.r.t. ${^θ}^{⋆}$ ;
$V$ and $- Δ V$ are both globally positive definite w.r.t. ${^θ}^{⋆}$ and $V$ is radially unbounded;
There exist class- $K$ functions $α_{1}, α_{2}, α_{3}$ such that $α_{2} (s) \leq α_{1} (ρ s) + α_{3} (s)$ holds near $s = 0$ , for some $ρ \in (0, 1)$ , and, for every $^θ$ near ${^θ}^{⋆}$ :

$α_{1} (∥^θ - {^θ}^{⋆} ∥) \leq V (^θ)$ $\leq α_{2} (∥^θ - {^θ}^{⋆} ∥)$ (6a)

$Δ V (^θ)$ $\leq - α_{3} (∥^θ - {^θ}^{⋆} ∥) .$ (6b)

Then, ${^θ}^{⋆}$ is

[noitemsep,topsep=0pt]
stable if 1 and 2 hold;
asymptotically stable if 1 and 3 hold;
globally asymptotically stable if 4 holds;
exponentially stable if 5 holds.

See supplementary section for proof. Notice that $Δ V$ was designed so that $Δ V_{k} = V_{k + 1} - V_{k}$ , where $Δ V_{k} ≜ Δ V ({^θ}_{k})$ and $V_{k} ≜ V ({^θ}_{k})$ . Furthermore, it is intended to represent a discrete-time equivalent to $˙ V (^θ (t)) = \frac{d}{d t} V (^θ (t))$ . The negative semidefiniteness simply translates to non-strict monotonicity of ${V_{k}}_{k \in Z_{+}}$ , which can be interpreted as a surrogate to ${ℓ(^θk)}k∈Z+$ .

4 Main Results: Convergence of EM

We now provide conditions that establish the different notions of Lyapunov stability explored thus far, for the dynamical system representation of the EM algorithm, and make appropriate conclusions in terms of the convergence of EM. To proceed, we first propose the natural candidate for a Lyapunov function in optimization, the function

V (^θ) ≜ ℓ (^θ) - ℓ ({^θ}^{⋆}) = log p ({^θ}^{⋆} | y) - log p (^θ | y),

(7)

where ${^θ}^{⋆}$ is a particular strict local maximum of interest (fixed for the remaining of this section), i.e. a particular point in the parameter space (state space) for which we seek convergence of EM. Since employing the Lyapunov theorem requires $V$ to be continuous, we make a mild assumption on the continuity of the incomplete-data posterior and then establish (non-asymptotic) stability.

Assumption 4.

$θ \mapsto log p (θ | y)$ is continuous.

Proposition 3.

Under Assumptions 1–4, any strict local maximum ${^θ}^{⋆}$ of the incomplete-data log-posterior $θ \mapsto log p (θ | y)$ is a stable equilibrium of the dynamical system (5) that represents EM.

Proof.

Consider the candidate Lyapunov function (7), defined for $^θ$ near ${^θ}^{⋆}$ . Clearly, $V$ is positive definite w.r.t. ${^θ}^{⋆}$ . Furthermore, notice that, since $Q (F (^θ),^θ) = max θ Q (θ,^θ)$ , then $Q (F (^θ),^θ) \geq Q (^θ,^θ)$ . Therefore, we have

log p

(F (^θ) | y) - D_{K L} [p_{_{X}} (\cdot | y,^θ) ∥ p_{_{X}} (\cdot | y, F (^θ))]      = d (F (^θ),^θ) \geq log p (^θ | y) - D_{K L} [p_{_{X}} (\cdot | y,^θ) ∥ p_{_{X}} (\cdot | y,^θ)]      = d (^θ,^θ) = 0,

(8)

where $d (θ,^θ)$ is defined in (4). In particular, we have


	$Δ V (^θ) = V (F (^θ)) - V (^θ) = log p (^θ \| y) - log p (F (^θ) \| y) \leq - d (F (^θ),^θ) \leq 0,$		(9a)

thus the result follows by invoking Lyapunov’s theorem. ∎

To establish asymptotic stability, it suffices that the inequality in (8) be strict, since through the same argument used for non-asymptotic stability, we thus have $Δ V (^θ) < 0$ for $∥^θ - {^θ}^{⋆} ∥ > 0$ sufficiently small. In order to achieve this, we make the following assumption.

Assumption 5.

$d (θ,^θ) = 0$ if and only if $θ =^θ$ , for $^θ$ near ${^θ}^{⋆}$ and arbitrary $θ$ .

This is also a relatively mild condition when ${^θ}^{⋆}$ is a “good” local minimizer of $ℓ (θ)$ . It follows, for instance, from a strong form of identifiability of the parameterized posterior latent distribution. It suffices that, for $^θ$ near ${^θ}^{⋆}$ , whenever $θ =^θ$ we have

P_{x \sim p (x | y,^θ)} [p (x | y, θ) = p (x | y,^θ)] = 1.

(10)

Failure to have Assumption 5 be satisfied near a particular strict local maximizer ${^θ}^{⋆}$ could result in that point not being being a fixed point of EM (i.e. equilibrium of EM’s dynamical system representation), let alone an asymptotically stable point and thus not locally convergent or even a limit point. Unfortunately, it is a well-documented behavior of the EM algorithm that its convergence properties and overall performance heavily rely on whether the initialization occurred near a “good” local maximizer of the incomplete-data log-likelihood or log-posterior [Figueiredo2004]. In the context of GMM clustering, for instance, Assumption 5 simply means that, after having sampled the GMM at hand, then the different posterior class probabilities will strictly depend on any perturbation to the parameters in the model (namely, the prior class probabilities and their corresponding means and covariance matrices).

4.1 Local and Global Convergence of EM

We now formally state and prove the local convergence of EM as a consequence of asymptotic Lyapunov stability. The proof is the straightforward culmination of previous section’s discussion.

Theorem 1 (Local Convergence).

Let ${^θ}^{⋆}$ a strict local maximizer of the incomplete-data posterior $θ \mapsto p (θ | y)$ and an isolated fixed point of EM. If Assumptions 1–5 hold, then ${^θ}^{⋆}$ is an asymptotically stable equilibrium of the dynamical system (5) that represents the EM algorithm (1), and thus EM is locally convergent to ${^θ}^{⋆}$ .

Proof.

Consider, once again, the candidate Lyapunov function (7), defined for $^θ$ near ${^θ}^{⋆}$ . Following the same strategy used in the proof of Proposition 3, we see that that $V$ is continuous and positive definite w.r.t. ${^θ}^{⋆}$ , but now $Δ V$ is not only negative semidefinite w.r.t. ${^θ}^{⋆}$ , but also negative definite, due to Assumption 5. ∎

In practice, since fixed points of EM must be stationary points of the log-posterior [Figueiredo2004], then a sufficient condition for this assumption to hold would be the continuous differentiability of $ℓ (^θ)$ for $^θ$ near ${^θ}^{⋆}$ , and that ${^θ}^{⋆}$ were an isolated stationary point. We now state the conditions that lead to global convergence of EM.

Theorem 2.

Suppose that ${θ : p (θ | y) > 0} = R^{p}$ and that the conditions of Theorem 1 hold, with Assumption 5 holding globally (i.e. $d (θ,^θ) = 0$ if and only if $θ =^θ$ for every $θ,^θ \in R^{p}$ ). If $ℓ (^θ) ≜ - log p (θ | y)$ is radially unbounded and ${^θ}^{⋆}$ is the only fixed point of EM, then ${^θ}^{⋆}$ is a globally asymptotically stable equilibrium of the dynamical system (5) that represents the EM algorithm (1), and thus EM is globally convergent to ${^θ}^{⋆}$ .

Proof.

We can reuse the argument in the proofs of Proposition 3 and Theorem 1 to check that condition 4 of the Lyapunov theorem holds with the additional assumptions. ∎

In particular, the unicity of fixed points follows, for the case of continuously differentiable incomplete-data log-posterior, for unimodal distributions. Furthermore, we require the support of the posterior to be the entire Euclidean space of appropriate dimension, and to vanish radially, (i.e. $p (θ | y) \to 0$ as $∥ θ ∥ \to \infty$ ). These last two conditions are largely technical and can be roughly circumvented in practice (for absolutely continuous posterior distributions). On the other hand, the unimodality rarely occurs and, in fact, EM can easily converge to saddle points or diverge.

4.2 Linear and Quadratic Convergence of EM

Returning to the assumptions that lead to exponential stability, and thus local convergence of EM, we now explore alternatives that can lead to linear and quadratic convergence of EM. We achieve this by further strengthening the strong identifiability Assumption 5. In addition, we want condition 5 of the Lyapunov theorem to be satisfied for the candidate Lyapunov function (7).

Assumption 6.

There exist some class- $K$ functions $α_{1}, α_{2}, α_{3}$ such that $α_{2} (s) \leq α_{1} (ρ s) + α_{3} (s)$ holds near $s = 0$ , for some $ρ \in (0, 1)$ , and for every $θ,^θ$ with ${^θ}^{⋆}$ near ${^θ}^{⋆}$

e^{- α_{2} (∥^θ - {^θ}^{⋆} ∥)} \leq \frac{p (^θ | y)}{p ({^θ}^{⋆} | y)} \leq e^{- α_{1} (∥^θ - {^θ}^{⋆} ∥)}, and, d (θ,^θ) \geq α_{3} (∥^θ - {^θ}^{⋆} ∥) .

(11)

This assumption sets the stage to state conditions for the linear convergence (in the sense of optimization) of the EM algorithm iterates.

Theorem 3 (Linear convergence of iterates).

Under the conditions of Theorem 1 and Assumption 6, ${^θ}^{⋆}$ is an exponentially stable equilibrium (with rate $ρ$ ) of the dynamical system (5) that represents the EM algorithm (1), and thus $∥ {^θ}_{k} - {^θ}^{⋆} ∥ \leq ∥ {^θ}_{0} - {^θ}^{⋆} ∥ \cdot ρ^{k} = O (ρ^{k})$ .

Similarly to what is noted in [Taylor2018], if $ℓ (^θ)$ is Lipschitz continuous, then we clearly have $ℓ ({^θ}_{k}) - ℓ ({^θ}^{⋆}) \leq L ∥ {^θ}_{0} -^θ^{⋆} ∥ \cdot ρ^{k} = O (ρ^{k})$ . In that sense, $p ({^θ}_{k} | y) \to p ({^θ}^{⋆} | y)$ actually converges quadratically, in the sense of optimization. Alternatively, it would have sufficed that $(α_{3} \circ α_{2}^{- 1}) (s) \geq (1 - μ) \cdot s$ for some $μ \in (0, 1)$ to achieve $ℓ ({^θ}_{k}) - ℓ ({^θ}^{⋆}) = O (μ^{k})$ . However, as discussed in Section 3.2, we can relax the condition in Lyapunov’s theorem required for exponential stability into something that leads to a notion of stability stronger than asymptotic stability but weaker than exponential stability, and which allows us to directly establish the Q-linear convergence of ${V_{k}}_{k \in Z_{+}}$ .

Theorem 4 (Quadratic convergence of posterior).

Under the conditions of Theorem 1, if there exists some $μ \in (0, 1)$ such that

p (^θ | y) \leq p ({^θ}^{⋆} | y) \cdot e^{- \frac{d (F (^θ),^θ)}{1 - μ}},

(12)

for every $^θ$ near ${^θ}^{⋆}$ , then $lim k \to \infty log p ({^θ}_{k} | y) = log p ({^θ}^{⋆} | y)$ with a Q-linear convergence rate upper bounded by $μ$ . In particular, we have $log p ({^θ}^{⋆} | y) - log p ({^θ}_{k} | y) \leq log (\frac{p ({^θ}^{⋆} | y)}{p ({^θ}_{0} | y)}) μ^{k} = O (μ^{k})$ .

Notice that (12) can be restated as $d (F (^θ),^θ) \geq (1 - μ) log (\frac{p ({^θ}^{⋆} | y)}{p (^θ | y)}),$ for $^θ$ near ${^θ}^{⋆}$ , to more closely resemble (11). Naturally, the main disadvantage of the last result is that checking the inequality (12) is likely to be virtually impossible in practice, given that it’s directly based on the EM dynamics $F$ . Furthermore, it has been widely observed that EM’s convergent rate is often sublinear, so the conditions in the previous results likely only hold in a few “good” local maximizers of the incomplete-data log-posterior. The last result focus on linear convergence of the posterior, by using $V (^θ) ≜ p ({^θ}^{⋆} | y) - p (^θ | y)$ as the new candidate Lyapunov function.

Theorem 5 (Linear convergence of the posterior).

Under the conditions of Theorem 1, if the concavity-like condition

p (F (^θ) | y) \geq μ p (^θ) + (1 - μ) p ({^θ}^{⋆} | y),

(13)

holds for every $^θ$ near ${^θ}^{⋆}$ , then $p ({^θ}^{⋆} | y) - p ({^θ}_{k} | y) \leq (p ({^θ}^{⋆} | y) - p ({^θ}_{0} | y)) \cdot μ^{k} = O (μ^{k})$ .

For the ease of readability, the proofs of Theorems 3-5 are detailed in the supplementary section.

4.3 Experimental validation of the convergence properties

We demonstrate the convergence properties of the MAP-EM algorithm on a general GMM with independent Gaussian priors on the unknown means. The details can be found in Appendix K, where we can see that, as the prior becomes more informative, convergence is achieved at faster rate. We also note that our convergence results readily apply to the MAP-EM algorithm over many distributions other than GMMs, and thus the wide range of applications it has been applied to.

5 Conclusion and Next Steps

In this paper, we addressed a gap in analyzing and designing optimization algorithms from the point of view of dynamical systems and control theory. Most recent literature largely focus on continuous-time

representations (via ordinary differential equations or inclusions) of general-purpose

gradient-based nonlinear optimization algorithms. However, we provide a unifying framework for the study of iterative optimization algorithms as discrete-time dynamical systems, and describe several relationships between forms of Lyapunov stability of state-space dynamical systems and convergence of optimization algorithms. In particular, we explored how exponential stability can be used to derive linear (or superlinear) convergence rates. We then narrowed this framework in detail to analyze convergence of the expectation-maximization (EM) algorithm for maximum a posteriori (MAP) estimation. Following first principles from dynamical systems stability theory, conditions for convergence of MAP-EM were developed, including conditions to show fast convergence (linear or quadratic), though EM often converges sublinearly.

The conditions we derived have a convenient statistical and information-theoretic interpretation and may thus be used in the future to design other novel EM-like algorithms with provable convergence rate. The conditions we derive would have been difficult to unveil without our approach, and thus we argue that a treatment similar to ours can we adopted for the convergence analysis of many other algorithms in machine learning. For future work, we believe that our approach can prove valuable in the design of EM-like algorithms for online estimation subject to a concept drift or data poisoning. In fact, by carefully designing a controller (input) on an otherwise unstable system, we can leverage a process known as stabilization to force different notions stability of dynamical systems that represent an optimization algorithm, and thus improve its convergence rate or robustness.

Acknowledgements

This work was supported by the Rensselaer-IBM AI Research Collaboration (http://airc.rpi.edu), part of the IBM AI Horizons Network (http://ibm.biz/AIHorizons).

References

Appendix

Appendix A Convergence via Subexponential Stability

As is often remarked when addressing in the context of EM, monotonicity is generally not enough to attain convergence. Nevertheless, the negative definiteness condition ensures asymptotic stability, and thus (local) convergence. Indeed, from an optimization perspective, the Lyapunov conditions for exponential stability are simply equivalent to strict monotonicity of ${V_{k}}_{k \in Z_{+}}$ , together with an absolute attainable lower bound of the surrogate $V (^θ)$ of $ℓ (^θ)$ at $^θ = {^θ}^{⋆}$ .

Note that from (6a) it follows that

Δ V (^θ) \leq - (α_{3} \circ α_{2}^{- 1}) (V (^θ))

(14)

for every $^θ$ near ${^θ}^{⋆}$ , which can be restated as

V_{k + 1} \leq (i d - α_{3} \circ α_{2}^{- 1}) (V_{k}),

(15)

for ${^θ}_{0}$ near ${^θ}^{⋆}$ , where $i d (s) ≜ s$ . Therefore, $V_{k} \leq (i d - α_{3} \circ α_{2}^{- 1})^{k} (V_{0})$ , where $α^{0} ≜ i d$ and $α^{k + 1} ≜ α \circ α^{k}$ . In particular, if $α_{i} (s) = a_{i} s^{p}$ with $a_{i} > 0$ ( $i = 1, 2, 3$ ) and $p > 0$ , then (14) and (15) become

Δ V_{k} \leq - (1 - μ) V_{k}

(16)

and

V_{k + 1} \leq μ V_{k},

(17)

respectively, where $μ ≜ 1 - a_{3} / a_{2} \in [0, 1)$ . In that case, we have $V_{k} \to 0$ as $k \to \infty$ , with a Q-linear convergence rate upper bounded by $μ$ . In particular, we have $V_{k} \leq V_{0} \cdot μ^{k} = O (μ^{k})$ . We further note that $α_{1}$ and $ρ$ do not directly influence the bound $μ$ . This observation is similar in spirit to Lemma 1 in [Aitken1994].

With these remarks into consideration, we note that if $V$ is continuous, positive definite w.r.t. ${^θ}^{⋆}$ , and (16) (equivalently, (17)) holds for $^θ$ near ${^θ}^{⋆}$ , then the linear rate $V_{k} = O (μ^{k})$ still holds, without necessarily having linear convergence of ${^θk}k∈Z+$ . Therefore, (14) without necessarily (6a) induces a weaker notion than exponential stability, but stronger than asymptotic stability. In fact, it coincides with the notion of $ℓ_{p}$ -stability for $α_{i} (s) = a_{i} s^{p^{i}}$ , with $p = p_{3} / p_{2} \geq 1$ [Lakshmikantham2002].

From an optimization perspective, we can leverage the ideas discussed in the previous paragraphs by noting that we are often concerned about the convergence rate in terms of the objective function or some meaningful surrogate of it, rather than directly the iterates in a numerical optimization scheme. This approach was implicit, for instance, in [Taylor2018, Vaquero2019].

Appendix B Relationships between stability and convergence

We represent the relationships between all the aforementioned concepts (local variants only, for the sake of simplicity) through the following diagram:

[arrows=Rightarrow] & exponentially stable [d][r] & linearly convergent [d]

& asymptotically stable [d][r] & locally convergent [d]

& stable point [d] & limit point [d]

& equilibrium [r, Leftrightarrow] & fixed point.

Appendix C Proof of Proposition 1

From Assumption 1, we have

p (y | x) = p (y | x, θ) = \frac{p (y | θ)}{p (x | θ)} p (x | y, θ) = \frac{1}{p (θ | x)} \frac{p (y)}{p (x)} p (θ | y) p (x | y, θ),

(18)

and therefore

p (θ | x) = \frac{p (y) / p (x)}{p (y | x)} p (θ | y) p (x | y, θ) = \frac{p (θ | y) p (x | y, θ)}{p (x | y)} .

(19)

Let $f_{1} (θ) \sim f_{2} (θ)$ indicate that $θ \mapsto f_{1} (θ) - f_{2} (θ)$ is constant. Then, we have

log p (θ | x) \sim log p (θ | y) + log p (x | y, θ) = log p (θ | y) - log (\frac{p (x | y,^θ)}{p (x | y, θ)})

(20)

for any fixed $^θ$ . Taking the expected value in $x \sim p (x | y,^θ)$ and attending to the definitions of EM (1) and KL divergence (Subsection 2.1), then the expression (2) follows. Finally, (3) readily follows by combining (1) and (2) and disregarding any additive terms that do not depend on $θ$ , since those will not affect the so-called M-step (maximization step) of the EM algorithm, i.e. the maximization in (1).

Appendix D Proof of Lemma 1

Let ${^θk}$ a convergent sequence, not necessarily generated by the EM algorithm (i.e. without necessarily having ${^θ}_{k + 1} = F ({^θ}_{k})$ ), with limit $lim k \to \infty {^θ}_{k} =^θ$ . Notice that


$Q (lim k \to \infty F ({^θ}_{k}),^θ)$	$= Q (lim k \to \infty F ({^θ}_{k}), lim k \to \infty {^θ}_{k})$	(21a)
	$= lim k \to \infty Q (F ({^θ}_{k}), {^θ}_{k})$	(21b)
	$\geq lim k \to \infty Q (F (^θ), {^θ}_{k})$	(21c)
	$= Q (F (^θ), lim k \to \infty {^θ}_{k})$	(21d)
	$= Q (F (^θ),^θ)$	(21e)
	$\geq Q (lim k \to \infty F ({^θ}_{k}),^θ),$	(21f)

where (21b) and (21d) both follow from the continuity of the $Q$ -function (Assumption 3). On the other hand, the inequalities (21c) and (21f) follow by noting that $Q (F ({^θ}_{k}), {^θ}_{k}) \geq F (θ, {^θ}_{k})$ and $Q (F (^θ),^θ) \geq Q (θ,^θ)$ for every $θ$ . In particular, they follow by choosing $θ = F (^θ)$ and $θ = lim k \to \infty F ({^θ}_{k})$ , respectively, and taking the limit $k \to \infty$ on both sides.

Therefore, we have

Q (lim k \to \infty F ({^θ}_{k}),^θ) = Q (F (^θ),^θ) = max θ Q (θ,^θ),

(22)

and thus $limk→∞F(^θk)∈argmaxθQ(θ,^θ)={F(^θ)}$ , which in turn makes $F$ continuous.

Appendix E Proof of Proposition 2

Let ${^θ}_{0}$ be such that $lim k \to \infty {^θ}_{k} = {^θ}^{⋆}$ , where ${^θk}k∈Z+$ was generated (1). Then,

F ({^θ}^{⋆}) = F (lim k \to \infty {^θ}_{k}) = lim k \to \infty F ({^θ}_{k}) = lim k \to \infty {^θ}_{k + 1} = {^θ}^{⋆},

(23)

where the second equality follows from the continuity of $F$ .

Appendix F Proof of Lemma 2

Let $ε^{(n)} > 0$ such that $ε^{(n)} \to 0$ as $n \to \infty$ and let $δ^{(n)} > 0$ be such that, if $∥ {^θ}_{0} - {^θ}^{⋆} ∥ \leq δ^{(n)}$ , then $∥ {^θ}_{k} - {^θ}^{⋆} ∥ \leq ε^{(n)}$ for every $k \in Z$ . Naturally, $0 < δ^{(n)} \leq ε^{(n)}$ , and thus $δ^{(n)} \to 0$ as $n \to \infty$ . Therefore, given any sequence ${^θ(n)0}n∈Z+$ such that $∥ {^θ}_{0}^{(n)} - {^θ}^{⋆} ∥ \leq δ^{(n)}$ , then ${^θ}_{0}^{(n)} \to {^θ}^{⋆}$ as $n \to \infty$ . Furthermore, $∥ F ({^θ}_{0}^{(n)}) - {^θ}^{⋆} ∥ = ∥ {^θ}_{1}^{(n)} - {^θ}^{⋆} ∥ \leq ε^{(n)}$ , and thus $F ({^θ}_{0}^{(n)}) \to {^θ}^{⋆}$ as $n \to \infty$ . From the continuity of $F$ , it follows that $F ({^θ}^{⋆}) = F (lim n \to \infty {^θ}_{0}^{(n)}) = lim n \to \infty F ({^θ}_{0}^{(n)}) = {^θ}^{⋆}$ .

Appendix G Proof of Lemma 3 (Lyapunov Theorem)

For the non-asymptotic and asymptotic stability part of the Lyapunov theorem, please refer to the proofs of Theorems 5.9.1–5.9.2 in [Agarwal2000], Corollary 4.8.1 and Theorem 4.8.3 in [Lakshmikantham2002], or Theorem 1.2 in [Bof2018]. For the global asymptotic stability, please refer to Theorem 5.9.8 in [Agarwal2000], Theorem 4.9.1. in [Lakshmikantham2002], or Theorem 1.4 in [Bof2018].

For the exponential stability, we can adapt the ideas in [Aitken1994] and the proof of Theorem 5.4 in [Bof2018]. More precisely, we have

α_{1} (∥ {^θ}_{k + 1} - {^θ}^{⋆} ∥) \leq V_{k + 1} = Δ V_{k} + V_{k} \leq α_{2} (∥ {^θ}_{k} - {^θ}^{⋆} ∥) - α_{3} (∥ {^θ}_{k} - {^θ}^{⋆} ∥) \leq α_{1} (ρ ∥ {^θ}_{k} - {^θ}^{⋆} ∥),

(24)

and thus $∥ {^θ}_{k + 1} - {^θ}^{⋆} ∥ \leq ρ ∥ {^θ}_{k} - {^θ}^{⋆} ∥$ . In other words, ${^θ}_{k} \to {^θ}^{⋆}$ as $k \to \infty$ Q-linearly, with a rate upper bounded by $ρ$ .

Appendix H Proof of Theorem 3 (Linear Convergence of Iterates)

Consider the candidate Lyapunov function $V (^θ) ≜ log p ({^θ}^{⋆} | y) - log p (^θ | y)$ , which is continuous due to Assumption 4. As seen in the proof of Proposition 3, we have $Δ V (^θ) \leq - d (F (^θ),^θ)$ . But then, from (11), we obtain (6b). Therefore, condition 5 of Lemma 3 is satisfied, and thus exponential stability is certified.

Appendix I Proof of Theorem 4 (Quadratic Convergence of Posterior)

Consider once again the continuous and positive definite candidate Lyapunov function $V (^θ) ≜ log p ({^θ}^{⋆} | y) - log p (^θ | y)$ . We will leverage the ideas in Appendix A. Indeed, notice that (12) can be rewritten as

p ({^θ}_{k} | y) \leq p ({^θ}^{⋆} | y) e^{- \frac{d (^θ_{k + 1}, {^θ}_{k})}{1 - μ}} .

(25)

Taking the log and rearranging terms, we get

d ({^θ}_{k + 1}, {^θ}_{k}) \leq (1 - μ) [log p ({^θ}^{⋆} | y) - log p ({^θ}_{k} | y)] = (1 - μ) V_{k} .

(26)

Finally, since $Δ V_{k} \leq - d ({^θ}_{k + 1}, {^θ}_{k}) \leq - (1 - μ) V)_{k}$ , then $Δ V_{k} \leq - (1 - μ) V_{k} ⟺ V_{k + 1} \leq μ V_{k}$ follows.

Appendix J Proof of Theorem 5 (Linear Convergence of Posterior)

This time, we consider the continuous and positive definite candidate Lyapunov function $V (^θ) ≜ p ({^θ}^{⋆} | y) - p (^θ | y)$ . Clearly, (13) can be rewritten as

p ({^θ}_{k + 1} | y) \geq μ p ({^θ}_{k} | y) + (1 - μ) p ({^θ}^{⋆} | y),

(27)

and thus

V_{k + 1} = p ({^θ}^{⋆} | y) - p ({^θ}_{k + 1} | y) \leq μ (p ({^θ}^{⋆} | y) - p ({^θ}_{k} | y) = μ V_{k} .

(28)

The result follows from the discussion in Appendix A.

Appendix K An Illustrative Example & Experiments

Let $y = {y^{(1)}, \dots, y^{(n)}}$ be some dataset that we wish to cluster. Let us assume that $y^{(i)} \in R^{d}$ is randomly, but independently, sampled from an unknown class $z^{(i)} \in {1, \dots, M}$ , which we thus seek to infer. For the sake of simplicity, consider that the class probabilities $P (z^{(i)} = m) = α_{m} \in (0, 1)$ such that


$α_{1} (∥^θ - {^θ}^{⋆} ∥) \leq V (^θ)$	$\leq α_{2} (∥^θ - {^θ}^{⋆} ∥)$	(6a)
$Δ V (^θ)$	$\leq - α_{3} (∥^θ - {^θ}^{⋆} ∥) .$	(6b)

A Dynamical Systems Approach for Convergence of the Bayesian EM Algorithm

Authors

Convergence of the Expectation-Maximization Algorithm Through Discrete-Time Lyapunov Stability Theory

On the Convergence of Bound Optimization Algorithms

Analysis of Gradient-Based Expectation-Maximization-Like Algorithms via Integral Quadratic Constraints

Optimization algorithms inspired by the geometry of dissipative systems

On Dissipative Symplectic Integration with Applications to Gradient-Based Optimization

A MAP approach for ℓ_q-norm regularized sparse parameter estimation using the EM algorithm

Optimization with Momentum: Dynamical, Control-Theoretic, and Symplectic Perspectives

This week in AI

1 Introduction

Contributions

2 Background: MAP-EM Algorithm

2.1 An Information-Theoretic Perspective

Assumption 1.

Proposition 1.

3 The Dynamical Systems Approach

Assumption 2.

3.1 Equilibrium in Dynamical Systems

Assumption 3.

Lemma 1.

Proposition 2.

3.2 Lyapunov Stability

Lemma 2.

Lemma 3 (Lyapunov theorem).

4 Main Results: Convergence of EM

Assumption 4.

Proposition 3.

Proof.

Assumption 5.

4.1 Local and Global Convergence of EM

Theorem 1 (Local Convergence).

Proof.

Theorem 2.

Proof.

4.2 Linear and Quadratic Convergence of EM

Assumption 6.

Theorem 3 (Linear convergence of iterates).

Theorem 4 (Quadratic convergence of posterior).

Theorem 5 (Linear convergence of the posterior).

4.3 Experimental validation of the convergence properties

5 Conclusion and Next Steps

Acknowledgements

References

Appendix

Appendix A Convergence via Subexponential Stability

Appendix B Relationships between stability and convergence

Appendix C Proof of Proposition 1

Appendix D Proof of Lemma 1

Appendix E Proof of Proposition 2

Appendix F Proof of Lemma 2

Appendix G Proof of Lemma 3 (Lyapunov Theorem)

Appendix H Proof of Theorem 3 (Linear Convergence of Iterates)

Appendix I Proof of Theorem 4 (Quadratic Convergence of Posterior)

Appendix J Proof of Theorem 5 (Linear Convergence of Posterior)

Appendix K An Illustrative Example & Experiments