A view of Estimation of Distribution Algorithms through the lens of Expectation-Maximization

TL;DR

This paper demonstrates that many Estimation of Distribution Algorithms can be interpreted as Monte Carlo Expectation-Maximization algorithms, providing a rigorous statistical framework for understanding and analyzing EDAs.

Contribution

It establishes a novel connection between EDAs and EM algorithms, offering a new theoretical perspective on their operation and analysis.

Findings

EDAs can be formulated as Monte Carlo EM algorithms

Exact EM limit corresponds to infinite samples in EDAs

Provides a rigorous statistical foundation for EDAs

Abstract

We show that a large class of Estimation of Distribution Algorithms, including, but not limited to, Covariance Matrix Adaption, can be written as a Monte Carlo Expectation-Maximization algorithm, and as exact EM in the limit of infinite samples. Because EM sits on a rigorous statistical foundation and has been thoroughly analyzed, this connection provides a new coherent framework with which to reason about EDAs.

Equations

$${\mathbf{\theta }}^{\ast }=\underset{\mathbf{\theta }}{\mathrm{argmax}}{\mathbb{E}}_{p(\mathbf{z}|\mathbf{\theta })}[f(\mathbf{z})],$$

$$\mbox{$\bm{\theta}^{(t+1)}=\operatorname*{argmax}_{\bm{\theta}}\sum_{i=1}^{N}W(f(\bm{z}_{i}))\log p(\bm{z}_{i}|\bm{\theta})$},$$

$$\mathcal{L}(\mathbf{\varphi })=\log p(\mathbf{x}|\mathbf{\varphi })=\log \int p(\mathbf{x},\mathbf{z}|\mathbf{\varphi })d\mathbf{z},$$

$$F(q,\mathbf{\varphi })={\mathbb{E}}_{q(\mathbf{z})}[\log p(\mathbf{x},\mathbf{z}|\mathbf{\varphi })]+H[q],$$

$$F_{MAP}(q,\mathbf{\varphi })\equiv {\mathbb{E}}_{q(\mathbf{z})}[\log p(\mathbf{x},\mathbf{z}|\mathbf{\varphi })+\log p_0(\mathbf{\varphi })]+H[q],$$

$$\hat{\mathbf{\theta }}$$

$$=\underset{\mathbf{\theta }}{\mathrm{argmax}}\log {\mathbb{E}}_{p(\mathbf{z}|\mathbf{\theta })}[f(\mathbf{z})]$$

$$=\underset{\mathbf{\theta }}{\mathrm{argmax}}{\mathcal{L}}_{EDA}(\mathbf{\theta }),$$

$${\mathcal{L}}_{EDA}(\mathbf{\theta })$$

$$=\log {\mathbb{E}}_{q(\mathbf{z})}\left[\frac{p(\mathbf{z}|\mathbf{\theta })f(\mathbf{z})}{q(\mathbf{z})}\right]$$

$$\ge {\mathbb{E}}_{q(\mathbf{z})}\left[\log (p(\mathbf{z}|\mathbf{\theta })f(\mathbf{z}))\right]+H[q]$$

$$=F(q,\mathbf{\theta }),$$

$$\tilde{p}(\mathbf{z}|\mathbf{\theta })=\frac{p(\mathbf{z}|\mathbf{\theta })f(\mathbf{z})}{{\int }_{\mathcal{Z}}p(\mathbf{z}|\mathbf{\theta })f(\mathbf{z})d\mathbf{z}},$$

$$q^{(t+1)}(\mathbf{z})=\frac{{\sum }_{i=1}^{N}f({\mathbf{z}}_i){\delta }_{{\mathbf{z}}_i}(\mathbf{z})}{{\sum }_{i=1}^{N}f({\mathbf{z}}_i)},$$

$${\mathbf{\theta }}^{(t+1)}=$$

$$=$$

$$p(\mathbf{z}|\mathbf{\theta })=h(\mathbf{z})\exp \left(\mathbf{\eta }(\mathbf{\theta }{)}^{T}T(\mathbf{z})-A(\mathbf{\theta })\right),$$

$${\mathbf{\theta }}^{(t+1)}=\frac{{\sum }_{i=1}^{N}f({\mathbf{z}}_i)T({\mathbf{z}}_i)}{{\sum }_{i=1}^{N}f({\mathbf{z}}_i)}.$$

$${\mathbf{\theta }}^{(t+1)}=(1-\gamma ){\mathbf{\theta }}^{(t)}+\gamma {\tilde{\mathbf{\theta }}}^{(t+1)},$$

$$p_0(\mathbf{\theta }|\mathbf{\lambda })=\exp \left({\mathbf{\lambda }}_1^T\mathbf{\eta }(\mathbf{\theta })-{\lambda }_{2}A(\mathbf{\theta })-B(\mathbf{\lambda })\right),$$

$${\hat{\mathcal{L}}}_{EDA}(\mathbf{\theta })=\log {\mathbb{E}}_{p(\mathbf{z}|\mathbf{\theta })}[f(\mathbf{z})p_0(\mathbf{\theta }|\mathbf{\lambda })],$$

$${\mathbf{\theta }}^{(t+1)}=\underset{\mathbf{\theta }}{\mathrm{argmax}}\sum _{i=1}^{N}f({\mathbf{z}}_i)\log p({\mathbf{z}}_i|\mathbf{\theta })p_0(\mathbf{\theta }|\mathbf{\lambda }).$$

$${\mathbf{\theta }}^{(t+1)}$$

$$=(1-\gamma ){\mathbf{\theta }}^{(t)}+\gamma {\tilde{\mathbf{\theta }}}^{(t+1)},$$

$${\nabla }_{\mathbf{\theta }}$$

$$={\mathbb{E}}_{p(\mathbf{z}|{\mathbf{\theta }}^{(t)})}[f(\mathbf{z}){\nabla }_{\mathbf{\theta }}\log p(\mathbf{z}|\mathbf{\theta }){|}_{\mathbf{\theta }={\mathbf{\theta }}^{(t)}}],$$

$${\mathbf{\theta }}^{(t+1)}={\mathbf{\theta }}^{(t)}+\alpha \sum _{i=1}^{N}f({\mathbf{z}}_i){\nabla }_{\theta }\log p({\mathbf{z}}_i|\mathbf{\theta }){|}_{\mathbf{\theta }={\mathbf{\theta }}^{(t)}},$$

$${\mathbf{\theta }}^{(t+1)}=$$

$${\mathbf{\theta }}^{(t+1)}=\underset{\mathbf{\theta }}{\mathrm{argmax}}\left[{\mathcal{L}}_{EDA}(\mathbf{\theta })-D_{KL}(\tilde{p}(\mathbf{z}|{\mathbf{\theta }}^{(t)})||\tilde{p}(\mathbf{z}|\mathbf{\theta }))\right],$$

$${\mathbf{\theta }}^{(t+1)}$$

$$\displaystyle+\langle\nabla_{\bm{\theta}}\mathcal{L}_{EDA}(\bm{\theta})\big{|}_{\bm{\theta}=\bm{\theta}^{(t)}},\bm{\theta}-\bm{\theta}^{(t)}\rangle$$

$$-D_{KL}(\tilde{p}(\mathbf{z}|{\mathbf{\theta }}^{(t)})||\tilde{p}(\mathbf{z}|\mathbf{\theta })).$$

$$D_{KL}$$

$$=\frac{1}{2}(\mathbf{\theta }-{\mathbf{\theta }}^{(t)}{)}^{\top }I({\mathbf{\theta }}^{(t)})(\mathbf{\theta }-{\mathbf{\theta }}^{(t)})+O\left((\mathbf{\theta }-{\mathbf{\theta }}^{(t)}{)}^3\right),$$

$${\mathbf{\theta }}^{(t+1)}$$

Full text

††copyright: rightsretained††conference: ; ;

1. Introduction

††

This is an expanded version of the paper by the same title published in the proceedings of the The Genetic and Evolutionary Computation Conference (GECCO), 2020

Estimation of Distribution Algorithms (EDAs) are a widely used class of algorithms designed to solve optimization problems of the form $\bm{z}^{*}=\operatorname*{argmax}_{\bm{z}\in\mathcal{Z}}f(\bm{z})$, where $f:\mathcal{Z}\rightarrow\mathbb{R}$ is a function over a space of discrete or continuous inputs, $\mathcal{Z}$. Instead of solving this objective directly, EDAs solve the related objective

[TABLE]

where $p(\bm{z}|\bm{\theta})$ is a probability density over $\mathcal{Z}$, parameterized by $D$ parameters $\bm{\theta}\in\mathbb{R}^{D}$. When $p(\bm{z}|\bm{\theta})$ has the capacity to represent point masses on the maxima of $f(\bm{z})$, then these two formulations have the same optimal values. Reasons for using the latter formulation of Equation 1 include convenience for derivative-free optimization (Hansen and Ostermeier, 2001), enabling of rigorous analysis of associated optimization algorithms(Zlochin et al., 2004), and the leveraging of a probabilistic formulation to incorporate auxiliary information (Brookes et al., 2019). In many cases this objective is further modified by applying a monotonic shaping function, $W(\cdot)$ to $f(\bm{z})$ (Wierstra et al., 2008, 2014; Yi et al., 2009), which alters convergence properties of algorithms to solve it, but does not change the optima.

The general algorithmic template for EDAs is as follows. Beginning with an initial parameter setting, $\bm{\theta}^{(0)}$, of the parametrized density, $p(\bm{z}|\bm{\theta})$, each iteration $t\in\{0,1,...T\}$ of an EDA generally consists of three steps:

1. (1)

Draw $N$ samples, $\{\bm{z}_{i}\}_{i=1}^{N}$, from $p(\bm{z}|\bm{\theta}^{(t)})$. 2. (2)

Evaluate $W(f(\bm{z}_{i}))$ for each $\bm{z}_{i}$. 3. (3)

Find a $\bm{\theta}^{(t+1)}$ that uses the weighted samples and corresponding function evaluations to move $p(\bm{z}|\bm{\theta})$ towards regions of $\mathcal{Z}$ that have large function values.

The last step is generally accomplished by attempting to solve

[TABLE]

which can be seen as a weighted maximum likelihood problem with weights $W(f(\bm{z}_{i}))$. We refer to this set of steps as a “core” EDA because it ties together most EDAs.

In the case of Covariance Matrix Adaptation (CMA-ES), this core algorithm is often modified in a variety of ways to improve performance. For example, samples from previous iterations may be used, directly or indirectly, in the last step, resulting in smoothing of the parameter estimates across iterations. Adaptive setting of parameter-specific step sizes, and “path evolution” heuristics (Hansen and Ostermeier, 2001; Ollivier et al., 2017; Brookes et al., 2019) are also common. These layers on top of the core EDA have generally been derived in a manner specific only to CMA-ES, and are not readily generalizable to other EDAs. For this reason, we restrict ourselves to the core EDA algorithm just described. However, as we shall see, our formalism also allows for a large class of EDA parameter smoothing variations which encompass many, but not all, variations of CMA-ES.

EDAs are also distinguished by the choice of parametrized density, $p(\bm{z}|\bm{\theta})$, which we refer to here as the “search model”.111This object is often simply referred to as the ”probability distribution” (Kern et al., 2004), however; we use ”search model” to distinguish it from other probability distributions used herein. Many EDAs use exponential family models for the search model, most commonly the multivariate Gaussian as in CMA-ES. However,

Bayesian networks (Ocenasek and Schwarz, 2012), Boltzmann machines (Hu and Hu, 2011; Shim et al., 2013), Markov Random Fields (Shakya and McCall, 2007), Variational Autoencoder (Brookes et al., 2019), and others (Kern et al., 2004) have also been used. Unless otherwise specified, our analysis in the following is quite general and makes no assumptions on the choice of search model.

Our contributions We show that a large class of EDAs—including, but not limited to, variants of CMA-ES—can be written as Monte Carlo (MC) Expectation-Maximization (EM) (Beal, 1998), and in the limit of infinite samples, as exact EM (Dempster et al., 1977). Because EM sits on a rigorous statistical foundation and has been thoroughly analyzed, this connection provides a new framework with which to reason about EDAs. Within this framework, we also show how many parameter-smoothed variants of CMA-ES can be written as Monte Carlo maximum a posteriori-EM (MAP-EM), which suggests a possible avenue to develop similar smoothing algorithms for EDAs with other search models. Finally, leveraging our EDA-EM connection, we provide a new perspective on how EDAs can be seen as performing approximate natural gradient descent.

1.1. Related work

The Information Geometry Optimization (IGO) framework (Ollivier et al., 2017) unifies a large class of approaches—including EDAs such as CEM (Rubinstein and Davidson, 1999) and CMA-ES (Hansen and Ostermeier, 2001)—for solving the objective in Equation 1 by discretizing a natural gradient flow. Instantiations of IGO are most readily seen as tantamount to using the score function estimator (sometimes referred to as the “log derivative trick”) on Equation 1, combined with natural gradient, as in Natural Evolution Strategies (Wierstra et al., 2008, 2014). The IGO framework does not connect to EM; therefore, IGO provides a complementary viewpoint to that presented herein.

Akimoto et al. (2012) show how CMA-ES with rank-$\mu$ updates, but without global step size and evolution paths, can be viewed as a natural evolution strategy. In this context, they briefly remark on how a simplified CMA-ES can be viewed as performing generalized EM with a partial-M step using a natural gradient. The technical underpinnings of this result are restricted specifically to CMA-ES because of the dependence on a Gaussian search model, and do not readily generalize to the broader EDA setting considered herein.

Staines et al. introduce the notion of variational optimization, by explicitly considering the fact that the optimum of the objective function in Equation 1 is a lower bound on the optimum of $f$. They also clearly delineate conditions under which the bound can be satiated, and the objective function is convex and has derivatives that exist (Staines and Barber, 2013); no connections to EDAs are made.

2. Background: Free Energy view of EM

††

This is an expanded version of the paper by the same title published in the proceedings of the The Genetic and Evolutionary Computation Conference (GECCO), 2020

Before deriving the connection between EM and EDAs, we first review some needed background on EM. EM is a well-known method for performing maximum likelihood parameter estimation in latent variable models that exploits the structure of the joint likelihood between latent and observed variables (Dempster et al., 1977). Intuitively, each E-step imputes the latent variables, and the subsequent M-step then uses these “filled in” data to do standard maximum likelihood estimation, typically in closed form. EM iterates between these E- and M- steps until convergence. We use the Free Energy interpretation of EM and its accompanying generalizations (Neal and Hinton, 1999) in order to reveal the connection between EDAs and EM.

Let $\bm{x}$ and $\bm{z}$ be observed and latent variables, respectively. The task of maximum likelihood estimation in a latent variable model is to find $\hat{\bm{\phi}}=\operatorname*{argmax}_{\bm{\phi}}\mathcal{L}(\bm{\phi})$ where

[TABLE]

for some model density $p$ parameterized by $\bm{\phi}$. In Neal and Hinton (1999), the authors define a function known as the free energy, given by

[TABLE]

where $q(\bm{z})$ is any probability density over the latent variables, $H[q]$ is the entropy of $q$, and the term preceding the entropy is known as the expected complete log-likelihood. The free energy lower bounds the log-likelihood, $\mathcal{L}(\bm{\phi})\geq F(q,\bm{\phi})$, and this bound is satiated only when $q(\bm{z})$ is equal to the true posterior, $q(\bm{z})=p(\bm{z}|\bm{x},\bm{\phi})$. If the true posterior cannot be calculated exactly, one may approximate it in one of several ways, two of which are described next.

In the first, the posterior is approximated by restricting it to a parameterized family of distributions, $q(\bm{z}|\bm{\psi})$, where both the parameters of the likelihood and the variational family must now be estimated. This is known as variational EM. Unless the true posterior lies in the parameterized class, which is not typically the case, the bound $\mathcal{L}(\bm{\phi})\geq F(q,\bm{\phi})$ cannot be satiated, leading to a variational gap given by $D_{KL}(q(\bm{z})||p(\bm{z}|\bm{x},\bm{\phi}))$, where $D_{KL}$ denotes the KL divergence.

In another approximation, one draws samples from the true posterior, $\bm{z}_{i}\sim p(\bm{z}|\bm{x},\bm{\phi}^{(t)})$, to estimate the expected complete log likelihood. This is known as Monte Carlo EM (MC-EM). This also induces a “gap” between the true and approximate posterior (because of the use of finitely many samples), which can be similarly computed as $D_{KL}(q(\bm{z})||p(\bm{z}|\bm{x},\bm{\phi}))$, and hence we will refer to it also as a variational gap. One major distinction between variational and MC-EM is that variational EM explicitly attempts to minimize the variational gap at each iteration through optimization, while MC-EM only implicitly closes the gap as the number of samples goes to infinity.

All of EM, MC-EM, and variational EM can be viewed as alternating coordinate descent on the free energy function  (Neal and Hinton, 1999). In particular, the alternating updates at iteration $t$ are given by:

• •

E-step: $q^{(t+1)}\leftarrow\operatorname*{argmax}_{q}F(q,\bm{\phi}^{(t)})$, that is, compute or estimate the posterior over the hidden variables. It can be shown that this is equivalent to minimizing the variational gap by solving $q^{(t+1)}(\bm{z})=\operatorname*{argmin}_{q}D_{KL}(q(\bm{z})||p(\bm{z}|\bm{x},\bm{\phi}^{(t)})$.

• •

M-step: $\bm{\phi}^{(t+1)}\leftarrow\operatorname*{argmax}_{\bm{\phi}}F(q^{(t+1)},\bm{\phi})$, which is equivalent to maximizing the expected complete log likelihood with respect to $\bm{\phi}$, $\bm{\phi}^{(t+1)}\leftarrow\operatorname*{argmax}_{\bm{\phi}}\mathbb{E}_{q^{(t+1)}(\bm{z})}[\log p(\bm{x},\bm{z}|\bm{\phi})]$.

When the variational gap can be driven to zero, as in exact EM, this procedure is guaranteed to never decrease the likelihood in Equation 3. Moreover, the convergence properties of exact EM to both local and global minima have been carefully studied (e. g., Dempster et al. (1977); Neal and Hinton (1999); Balakrishnan et al. (2017)). Rigorous generalizations of EM to partial E- and M-steps also emerge naturally from this viewpoint (Neal and Hinton, 1999).

One can also use the free energy viewpoint of EM to optimize a maximum a posteriori (MAP) objective given by

$\mathcal{L}_{MAP}(\bm{\phi})\equiv\log p(\bm{x}|\bm{\phi})+\log p_{0}(\bm{\phi})$, where $p_{0}$ is some prior distribution over parameters. This yields a corresponding free energy,

[TABLE]

upon which one can perform the same coordinate descent algorithm just described. This formulation is referred to as MAP-EM. We will draw connections between MAP-EM and EDAs with smoothed parameter updates such as those commonly used in CMA-ES.

3. Formal connection between EM and EDAs

We are now ready to use the EM framework presented in Section 2, to show that EDAs using the update rule in Equation 2 can be written as MC-EM, and as exact EM in the infinite sample limit. We then show that generalizations of Equation 2 that allow for parameter smoothing—such as CMA-ES’s smoothing with the previous estimate and use of an evolutionary path—can be readily incorporated into the EM-EDA connection by using the MAP-EM framework. Finally, we highlight a connection between EM and standard gradient-based optimization.

3.1. Derivation of EDAs as MC-EM

††

This is an expanded version of the paper by the same title published in the proceedings of the The Genetic and Evolutionary Computation Conference (GECCO), 2020

As described in the introduction, EDAs seek to solve the objective defined in Equation 1,

[TABLE]

where $f(\bm{z})$ is the black-box function to be optimized, $p(\bm{z}|\bm{\theta})$ is what we refer to as the search model, parameterized by $\bm{\theta}$, and we define $\mathcal{L}_{EDA}(\bm{\theta})\equiv\log\mathbb{E}_{p(\bm{z}|\bm{\theta})}[f(\bm{z})]$. This expression can be thought of as an EDA equivalent to the ‘log marginal likelihood’ in EM, only without any observed data, $\bm{x}$.

Some EDAs monotonically transform $f(\bm{z})$ with a shaping function, $W(\cdot)$, which may be, for example, an exponential (Peters and Schaal, 2007), a

quantile-based transformation (Ollivier et al., 2017; Rubinstein and Davidson, 1999), or a cumulative density function (CDF) (Brookes et al., 2019). Although this transformation does not change the optima, it may alter the optimization dynamics. Often this shaping function is changed at each iteration in a sample-dependent, adaptive manner (which links these methods to annealed versions of EM and VI (Hofmann, 2001; Mandt et al., 2016)). In such a setting, the connection that we will show between EDAs and EM holds within each full iteration. For notational simplicity, we drop the $W(\cdot)$ and assume that $f(\bm{z})$ has already been transformed. We additionally assume that the transformation is such that $f(\bm{z})\geq 0$ for all $\bm{z}\in\mathcal{Z}$.

To link Equation 8 to EM, we introduce a probability density, $q(\bm{z})$, which allows us to derive a lower bound on $\mathcal{L}_{EDA}$ using Jensen’s inequality:

[TABLE]

where $F$ is the same free energy function appearing in the previous section on EM, except that the complete likelihood is replaced with the term $f(\bm{z})p(\bm{z}|\bm{\theta})$. When $f(\bm{z})p(\bm{z}|\bm{\theta})$ is normalizable, then it can be shown that there is an EDA “variational gap” given by $F(q,\bm{\theta})-\mathcal{L}_{EDA}=-D_{KL}(q(\bm{z})||\tilde{p}(\bm{z}|\bm{\theta}))$, where we define the tilted density,

[TABLE]

which is the EDA counterpart to the exact posterior in EM. We can now construct a coordinate ascent algorithm on the free energy defined in Equation 12 that mirrors the EM algorithm. In particular, this algorithm iterates between E-steps that solve

$q^{(t+1)}\leftarrow\operatorname*{argmin}_{q}D_{KL}(q(\bm{z})||\tilde{p}(\bm{z}|\bm{\theta}^{(t)}))$, and M-steps that solve

$\theta^{(t+1)}\leftarrow\operatorname*{argmax}_{\theta}\mathbb{E}_{q^{(t+1)}(\bm{z})}\left[\log(p(\bm{z}|\bm{\theta})f(\bm{z}))\right]$. To make the precise connection between practically implemented EDA and EM, we introduce a particular approximate ‘posterior’ for the E-step that is given by a mixture of weighted particles:

[TABLE]

where $\{\bm{z}_{i}\}_{i=1}^{N}$ are samples drawn from $p(\bm{z}|\bm{\theta}^{(t)})$, as in EDAs. Using this posterior approximation, the M-step amounts to solving the objective:

[TABLE]

which is exactly the generic EDA update step defined in Equation 2.

From this exposition it also becomes clear that EDAs can be thought of as performing a type of MC-EM that uses Importance Sampling (MacKay, 2003) rather than direct sampling from the posterior. In particular, the EDA sampling procedure uses proposal distribution, $p(\bm{z}|\theta^{(t)})$, (sampled in the E-step), and importance weights

$p(\bm{z}|\bm{\theta})f(\bm{z})/p(\bm{z}|\bm{\theta})=f(\bm{z})$ to estimate the expected “complete log likelihood", $\mathbb{E}_{q^{(t+1)}(\bm{z})}\left[\log(p(\bm{z}|\bm{\theta})f(\bm{z}))\right]$, in the M-step.

This is our main result, as it shows that many EDAs can be viewed as an EM algorithm that uses the particle-based posterior approximation given by Equation 14. For any $\bm{z}$, we have

$q^{(t+1)}(\bm{z})\overset{p}{\to}\tilde{p}(\bm{z}|\bm{\theta}^{(t)})$ as $n\to\infty$ by the law of large numbers and Slutsky’s theorem (Grimmett and Stirzaker, 2001). In this limit of infinite particles, the approximate posterior matches the tilted distribution—the EDA equivalent to the “exact posterior”—and our algorithm inherits the same guarantees as exact EM, such as guaranteed improvement of the objective function at each iteration, as well as local and global convergence properties (Neal and Hinton, 1999; Balakrishnan et al., 2017; Salakhutdinov et al., 2003).

In many cases, EDAs use a member of an exponential family for the search model. Letting $\bm{\theta}$ denote the expectation parameters, then the search model has the density:

[TABLE]

where $h$ is the base measure, $\bm{\eta}(\bm{\theta})$ are the natural parameters, $T(\bm{z})$ are sufficient statistics and $A(\bm{\theta})$ is the log-partition function. Then the update of Equation 17 takes the simple form of a weighted maximum likelihood estimate for the exponential family:

[TABLE]

Next we show how exponential family search models can be used to connect parameter-smoothed EDAs to MAP-EM.

3.2. Smoothed EDAs as MAP-EM

In CMA-ES, the parameter updates are smoothed between iterations in a number of ways (Hansen et al., 2003). For example, the covariance estimate is typically smoothed with the previous estimate. Additionally, it may be further smoothed using a rank one covariance matrix obtained from the “evolution path” that the algorithm has recently taken.222See the two terms in Equation 11 in Hansen et al. (2003)). Next, we consider these types of smoothing in detail. However, any smoothing update, or combination thereof, can be similarly derived by adjusting the form of the prior distribution. A benefit of viewing CMA-ES updates in this general form is that it becomes more straightforward to determine how to do similar types of smoothing for EDAs without Gaussian search models.

When smoothing with the previous parameter estimate, the updates can be written as

[TABLE]

where $\tilde{\bm{\theta}}^{(t+1)}$ is the solution to Equation 17 and $\gamma$ is a hyperparameter that controls the amount of smoothing. In the case where the search model is a member of an exponential family, as is defined in Equation 18, we will show that the smoothed update of Equation 20 is equivalent to a particular MAP-EM update that uses the tilted density of Equation 13 as the approximate posterior. To see this, consider the conjugate prior to the exponential family,

[TABLE]

where $\bm{\lambda}=(\bm{\lambda}_{1},\lambda_{2})$, $B(\bm{\lambda})$ is the log-partition function of the prior, and we have assumed that the base measure is constant. We now consider the modified EDA objective:

[TABLE]

which is analogous to the MAP objective defined in Section 2. It can be shown that by replacing Equation 9 with this modified objective, performing analogous steps as those in Equations (10-12) and (15-17), and using the same definition of the tilted density and approximate posterior as in Equations 13 and 14, respectively, we arrive at the update equation:

[TABLE]

If we now let $\lambda_{2}=\nicefrac{{1}}{{\gamma}}-1$ and allow $\bm{\lambda}_{1}$ to change every iteration as $\bm{\lambda}_{1}=(\nicefrac{{1}}{{\gamma}}-1)\bm{\theta}^{(t)}$, it can then be shown that the solution to this objective is given by:

[TABLE]

where $\tilde{\bm{\theta}}^{(t+1)}=\frac{\sum_{i=1}^{N}f(\bm{z}_{i})T(\bm{z}_{i})}{\sum_{i=1}^{N}f(\bm{z}_{i})}$ is Equation 19, the solution of original EDA objective. Since the update in Equation 24 is identical to the smoothed update of Equation 20, we can see that smoothing can be viewed as a consequence of performing MAP-EM in the context of EDAs, with a particular choice of prior parameters. See the Appendix A.2 for a full derivation.

3.3. A continuum between EM-like algorithms and gradient descent

Instead of using EDAs, one could alternatively solve the objective in Equation 1 directly with gradient descent (or natural gradient descent) where updates are given by $\bm{\theta}^{(t+1)}=\bm{\theta}^{(t)}+\alpha\nabla\mathcal{L}_{EDA}(\bm{\theta}^{(t)})$, and $\alpha$ is a step size parameter. This is known as REINFORCE (Williams, 1992) in the RL community and is related to IGO-ML (Ollivier et al., 2017). Typically, one cannot compute the gradient term exactly, in which case the “log derivative trick”,

[TABLE]

is often combined with Monte Carlo estimation to arrive at an update:

[TABLE]

where $\{\bm{z}_{i}\}_{i=1}^{N}$ are samples drawn from $p(\bm{z}|\bm{\theta}^{(t)})$. Notably, this update does not require the gradient of $f(\bm{z})$ (Recht, 2018).

We can connect the gradient-based optimization in Equation 27 to the EDA “EM” optimization presented in Section 3.1. In particular, consider the partial M-step version of EM, where instead of fully solving for the new parameter, one instead performs one gradient step —so called “first-order” EM (Neal and Hinton, 1999; Balakrishnan et al., 2017). Instantiating first-order EM into EDA “EM” by partially solving the objective in Equation 17 with a single gradient step (with step size $\alpha$) will result in the same update as Equation 27. In other words, performing “first-order" EM in the EDA context is identical to performing gradient descent on Equation 1. As we can see, as one performs more gradient steps within an M-step, one traces out a continuum of methods between gradient descent on the original objective and EM.

4. A unifying view of EDAs as natural gradient descent

In the next sections, we provide an illustrative example of how the connection between EDAs and EM can be leveraged. In particular, we use the framework to provide a new perspective on how EDAs can be seen as approximating natural gradient descent.

Akimoto et al. (2010) showed that CMA-ES, a particular EDA, can be viewed as approximate natural gradient descent (NGD) . This result also emerges from the IGO framework of Ollivier et al. (2017), who additionally showed that CEM and several other EDAs can be viewed as NGD. Malagò et al. (2013) also show how EDAs using MRF for the search model can be interpreted as NGD. The relationship we derived between EDA and EM provides a unifying perspective on these connections, which also allows for them to be readily generalized to any EDA that can be written as EM. Next we develop this unifying view.

††

This is an expanded version of the paper by the same title published in the proceedings of the The Genetic and Evolutionary Computation Conference (GECCO), 2020 It does not appear widely known that EM can be viewed as approximate NGD. Thus we next synthesize a series of known, but independent, results to show this. As a method for maximizing the log-likelihood function, Chrétien and Hero showed that EM can be formulated as the proximal point method (PPM) with the reverse KL divergence (Chrétien and Hero III, 1998; Chrétien and Hero, 2000). Specifically, this formulation shows that in our context, each update

[TABLE]

is equivalent to

[TABLE]

which is a PPM (see also Appendix B). This equivalence requires the exact posterior (tilted distribution), $\tilde{p}$, and therefore does not technically hold for EDA in practice, which is more akin to MC-EM. Nevertheless, the connection exposes an interesting new perspective on EDAs in the asymptotic regime, as we shall see.

Next, consider the mirror descent update, which is a linearization of the PPM (Beck and Teboulle, 2003):

[TABLE]

Mirror descent with KL divergence is equivalent to NGD on the dual Reimannian manifold (Raskutti and Mukherjee, 2015). Even when the original and dual manifolds differ, the two algorithms can be related: using a second-order Taylor series approximation of the KL divergence yields an approximate mirror descent update that is precisely NGD. Specifically, replacing the KL divergence in mirror descent with the second-order approximation as follows,

[TABLE]

where $I$ is the Fisher information metric (Martens, 2014), yields the update

[TABLE]

which can be shown to be equivalent to the NGD update

[TABLE]

See Appendix C for further details.

In summary: EM is a PPM with KL divergence; mirror descent is a linearization of the PPM; and a second-order approximation of the KL divergence in mirror descent yields NGD. Thus in the infinite sample limit, EDAs can be seen as approximate NGD.

Note that Sato (2001), and later others, showed the related result that for exponential family models with hidden variables and a conjugate prior (to the joint observed and hidden model), classical (non-amortized) mean-field variational inference performs NGD (Hensman et al., 2012; Hoffman et al., 2013). In contrast, our result applies to any exact EM, but only connects it to approximate NGD.

An interesting side note is that Salimbeni et al. (2018) demonstrated empirically that the best step size for NGD, according to Brent line searches, increases over iterations to 1.0, for Gaussian approximations of several models. Their result aligns nicely with our results here that EM, viewed as approximate NGD, has a step size of precisely 1.0.

5. Discussion

We have shown a novel connection between a broad class of EDAs and EM. Additionally, we have presented an illustrative example of insight from EM that can be applied to EDAs by way of this connection. Specifically, we presented a new connection between NGD and EDAs, by way of EM.

In light of EDAs as EM, one can ask whether this viewpoint sheds any light on choices of the EDA search model that may be more or less optimal. Suppose for a moment that rather than a sample-based approximate posterior, one instead used a parametric posterior, as in standard VI. In such a case, if the search model were in some sense “conjugate" to $f(\bm{z})$, then the posterior could take on an analytical form. When might such an approach make sense? Suppose one were using EDAs to perform protein design (Brookes et al., 2019), that $f(\bm{z})$ represented protein stability, and that $f(\bm{z})$ could be approximated by using an exponential model of linear additive marginal and pairwise terms (Bryngelson and Wolynes, 1987). Then if one used a Potts model (Marks et al., 2011) as the search model, the exact posterior could also be written in this form, and a “parametric" EDA could be pursued.

Given the tight connection between EDAs with adaptive shape functions, and annealed versions of EM and VI (Hofmann, 2001; Mandt et al., 2016), we believe that these latter approaches could benefit from the simple and robust implicit annealing schemes found in the EDA literature, which arise not from an annealing schedule, but from simple and effective quantile transformations and the like.

Appendix A Derivations in the case of Exponential Family search models

In Equations 19 and 23 we present exact update equations for EDAs when the search model is a member of an exponential family. We derive these updates here.

A.1. EDA Update

For the in update, Equation 19, we use the definition of the exponential family search model density in Equation 18 in the EDA update objective in Equation 17:

[TABLE]

To solve this, we now take the gradient of the argument with respect to $\bm{\theta}$:

[TABLE]

where we have used the chain rule and the properties of the log-partition function (i.e. cumulant generating function) to equate $\nabla_{\bm{\theta}}A(\bm{\theta})=\nabla_{\bm{\theta}}\bm{\eta}(\bm{\theta})\mathbb{E}_{\bm{\eta}}[T(\bm{z})]$. Now, recognizing that $\bm{\theta}$ are expectation parameters, so $\bm{\theta}=\mathbb{E}_{\bm{\eta}}[T(\bm{z})]$, and setting Equation 42 equal to zero, we arrive at the update in Equation 19:

[TABLE]

A.2. MAP Update

††

This is an expanded version of the paper by the same title published in the proceedings of the The Genetic and Evolutionary Computation Conference (GECCO), 2020

For the update in Equation 23, we use the definition of the exponential family search model density in Equation 18, and the prior in Equation 21, in the MAP-EDA update objective, Equation 22:

[TABLE]

Taking the gradient of the argument with respect to $\bm{\theta}$ and setting it equal to zero:

[TABLE]

where again we use $\bm{\theta}=\mathbb{E}_{\bm{\eta}}[T(\bm{z})]$. Now, we can recognize $\tilde{\bm{\theta}}^{(t+1)}=\frac{\sum_{i=1}^{N}f(\bm{z}_{i})T(\bm{z}_{i})}{\sum_{i=1}^{N}f(\bm{z}_{i})}$ is the solution to the original EDA objective, and letting $\lambda_{2}=\nicefrac{{1}}{{\gamma}}-1$ and $\bm{\lambda}_{1}=(\nicefrac{{1}}{{\gamma}}-1)\bm{\theta}^{(t)}$ we arrive at:

[TABLE]

which is the desired smoothed update.

Appendix B Equivalence between Expectation-Maximization and the proximal point method

In Section 6, we discuss the equivalence between EM and using the proximal point method (PPM) with reverse KL divergence to maximize the log-likelihood function (Chrétien and Hero III, 1998; Chrétien and Hero, 2000). Here, we provide the derivation of this equivalence, beginning with a formulation of the PPM:

[TABLE]

Note that the last line is precisely the update rule for EM.

Appendix C Details of natural gradient descent

Recall from Section 6 that replacing the log-likelihood with its linearization and the KL-divergence with its second-order Taylor series approximation in the PPM formulation of EM yields the update

[TABLE]

One can further derive a closed-form update by taking the derivative of the $\operatorname*{argmax}$ argument on the right-hand side, setting it equal to zero as per first-order optimality conditions, and solving:

[TABLE]

This establishes the equivalence between the standard natural gradient descent update rule and the mirror descent update derived using the PPM formulation of EM as a starting point.

Acknowledgements.

We thank Ben Recht for helpful discussion, David Duvenaud, Roger Grosse and James Hensman for pointers to variational Bayes as NGD, Petros Koumoutsakos for encouraging us to think about how to incorporate smoothing, Nikolaus Hansen for pointing us to relevant work, Youhei Akimoto for replies to questions about his papers, and Meng Xiangming for catching technical typos. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.