t\hskip 0.50003pt$ equal to $1$. However, it is well-known, both theoretically
and experimentally [BLS07, Han06b, WAS04], that standard CEM loses diversity
too fast in many situations. The usual solution [dBKMR05] is to reduce the
learning rate (smoothed CEM, (25)), but this breaks the reparametrization
invariance.

On the other hand, the IGO flow can be seen as a _maximum likelihood update
with infinitesimal learning rate_ (Theorem 13). This interpretation allows to
define a particular IGO algorithm, the IGO-ML (Definition 14): it performs a
maximum likelihood update with an arbitrary learning rate, and keeps the
reparametrization invariance. It coincides with CEM when the learning rate is
set to $1$, but it differs from smoothed CEM by the exchange of the order of
argmax and averaging (compare (23) and (25)). We argue that this new algorithm
may be a better way to reduce the learning rate and achieve smoothing in CEM.

Standard CEM can lose diversity, yet is a particular case of an IGO algorithm:
this illustrates the fact that reasonable values of the learning rate $\hskip
0.50003pt\delta t\hskip 0.50003pt$ depend on the parametrization. We have
studied this phenomenon in detail for various Gaussian IGO algorithms (Section
4.2).

Why would a smaller learning rate perform better than a large one in an
optimization setting? It might seem more efficient to jump directly to the
maximum likelihood estimate of currently known good points, instead of
performing a slow gradient ascent towards this maximum.

However, optimization faces a “moving target”, contrary to a learning setting
in which the example distribution is often stationary. Currently known good
points are likely not to indicate the _position_ at which the optimum lies,
but, rather, the _direction_ in which the optimum is to be found. After an
update, the next elite sample points are going to be located somewhere new. So
the goal is certainly not to settle down around these currently known points,
as a maximum likelihood update does: by design, CEM only tries to reflect
status-quo (even for $N=\infty$), whereas IGO tries to move somewhere. When
the target moves over time, a progressive gradient ascent is more reasonable
than an immediate jump to a temporary optimum, and realizes a kind of time
smoothing.

This phenomenon is most clear when the number of sample points is small. Then,
a full maximum likelihood update risks losing a lot of diversity; it may even
produce a degenerate distribution if the number of sample points is smaller
than the number of parameters of the distribution. On the other hand, for
smaller $\hskip 0.50003pt\delta t\hskip 0.50003pt$, the IGO algorithms do, by
design, try to maintain diversity by moving as little as possible from the
current distribution $P_{\theta}$ in Kullback–Leibler divergence. A full ML
update disregards the current distribution and tries to move as close as
possible to the elite sample in Kullback–Leibler divergence [dBKMR05], thus
realizing maximal diversity loss. This makes sense in a non-iterated scenario
but is unsuited for optimization.

##### Diversity and multiple optima.

The IGO framework emphasizes the relation of natural gradient and diversity:
we argued that IGO provides minimal diversity change for a given objective
function increment. In particular, provided the initial diversity is large,
diversity is kept at a maximum. This theoretical relationship has been
confirmed experimentally for restricted Boltzmann machines.

On the other hand, using the vanilla gradient does not lead to a balanced
distribution between the two optima in our experiments. Using the vanilla
gradient introduces hidden arbitrary choices between those points (more
exactly between moves in $\Theta$-space). This results in loss of diversity,
and might also be detrimental at later stages in the optimization. This may
reflect the fact that the Euclidean metric on the space of parameters,
implicitly used in the vanilla gradient, becomes less and less meaningful for
gradient descent on complex distributions.

IGO and the natural gradient are certainly relevant to the well-known problem
of exploration-exploitation balance: as we have seen, arguably the natural
gradient realizes the best increment of the objective function with the least
possible change of diversity in the population.

More generally, IGO tries to learn a model of where the good points are,
representing all good points seen so far rather than focusing only on some
good points; this is typical of machine learning, one of the contexts for
which the natural gradient was studied. The conceptual relationship of IGO and
IGO-like optimization algorithms with machine learning is still to be explored
and exploited.

## Summary and conclusion

We sum up:

  * •

The information-geometric optimization (IGO) framework derives from invariance
principles and allows to build optimization algorithms from any family of
distributions on any search space. In some instances (Gaussian distributions
on ${\mathbb{R}}^{d}$ or Bernoulli distributions on $\\{0,1\\}^{d}$) it
recovers versions of known algorithms (CMA-ES or PBIL); in other instances
(restricted Boltzmann machine distributions) it produces new, hopefully
efficient optimization algorithms.

  * •

The use of a quantile-based, time-dependent transform of the objective
function (Equation (3)) provides a rigorous derivation of rank-based update
rules currently used in optimization algorithms. Theorem 4 uniquely identifies
the infinite-population limit of these update rules.

  * •

The IGO flow is singled out by its equivalent description as an infinitesimal
weighted maximal log-likelihood update (Theorem 13). In a particular
parametrization and with a step size of $1$, it recovers the cross-entropy
method (Corollary 16).

  * •

Theoretical arguments suggest that the IGO flow minimizes the change of
diversity in the course of optimization. In particular, starting with high
diversity and using multimodal distributions may allow simultaneous
exploration of multiple optima of the objective function. Preliminary
experiments with restricted Boltzmann machines confirm this effect in a simple
situation.

Thus, the IGO framework is an attempt to provide sound theoretical foundations
to optimization algorithms based on probability distributions. In particular,
this viewpoint helps to bridge the gap between continuous and discrete
optimization.

The invariance properties, which reduce the number of arbitrary choices,
together with the relationship between natural gradient and diversity, may
contribute to a theoretical explanation of the good practical performance of
those currently used algorithms, such as CMA-ES, which can be interpreted as
instantiations of IGO.

We hope that invariance properties will acquire in computer science the
importance they have in mathematics, where intrinsic thinking is the first
step for abstract linear algebra or differential geometry, and in modern
physics, where the notions of invariance w.r.t. the coordinate system and so-
called gauge invariance play central roles.

## Acknowledgements

The authors would like to thank Michèle Sebag for the acronym and for helpful
comments. Y. O. would like to thank Cédric Villani and Bruno Sévennec for
helpful discussions on the Fisher metric. A. A. and N. H. would like to
acknowledge the Dagstuhl Seminar No 10361 on the Theory of Evolutionary
Computation (http://www.dagstuhl.de/10361) for inspiring their work on natural
gradients and beyond. This work was partially supported by the
ANR-2010-COSI-002 grant (SIMINOLE) of the French National Research Agency.

## Appendix: Proofs

### Proof of Theorem 4 (Convergence of empirical means and quantiles)

Let us give a more precise statement including the necessary regularity
conditions.

######  Proposition 18.

Let $\theta\in\Theta$. Assume that the derivative $\frac{\partial\ln
P_{\theta}(x)}{\partial\theta}$ exists for $P_{\theta}$-almost all $x\in X$
and that $\mathbb{E}_{P_{\theta}}\left|\mskip 1.0mu\frac{\partial\ln
P_{\theta}(x)}{\partial\theta}\right|^{2}<+\infty$. Assume that the function
$w$ is non-decreasing and bounded.

Let $(x_{i})_{i\in{\mathbb{N}}}$ be a sequence of independent samples of
$P_{\theta}$. Then with probability $1$, as $N\to\infty$ we have

$\frac{1}{N}\sum_{i=1}^{N}\widehat{W}^{f}(x_{i})\frac{\partial\ln
P_{\theta}(x_{i})}{\partial\theta}\to\int W_{\theta}^{f}(x)\,\frac{\partial\ln
P_{\theta}(x)}{\partial\theta}\,P_{\theta}({\mathrm{d}}x)$

where

$\widehat{W}^{f}(x_{i})=w\left(\frac{\mathrm{rk}_{N}(x_{i})+1/2}{N}\right)$

with $\mathrm{rk}_{N}(x_{i})=\\#\\{1\leqslant j\leqslant
N,f(x_{j})<f(x_{i})\\}$. (When there are $f$-ties in the sample,
$W^{f}(x_{i})$ is defined as the average of $w((r+1/2)/N)$ over the possible
rankings $r$ of $x_{i}$.)

###### Proof.

Let $g:X\to{\mathbb{R}}$ be any function with
$\mathbb{E}_{P_{\theta}}g^{2}<\infty$. We will show that
$\frac{1}{N}\sum\widehat{W}^{f}(x_{i})g(x_{i})\to\int
W_{\theta}^{f}(x)g(x)\,P_{\theta}({\mathrm{d}}x)$. Applying this with $g$
equal to the components of $\frac{\partial\ln P_{\theta}(x)}{\partial\theta}$
will yield the result.

Let us decompose

$\displaystyle\frac{1}{N}\sum\widehat{W}^{f}(x_{i})g(x_{i})=\frac{1}{N}\sum
W_{\theta}^{f}(x_{i})g(x_{i})+\frac{1}{N}\sum(\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i}))g(x_{i}).$

Each summand in the first term involves only one sample $x_{i}$ (contrary to
$\widehat{W}^{f}(x_{i})$ which depends on the whole sample). So by the strong
law of large numbers, almost surely $\frac{1}{N}\sum
W_{\theta}^{f}(x_{i})g(x_{i})$ converges to $\int
W_{\theta}^{f}(x)g(x)\,P_{\theta}({\mathrm{d}}x)$. So we have to show that the
second term converges to $0$ almost surely.

By the Cauchy–Schwarz inequality, we have

$\left|\mskip
1.0mu\frac{1}{N}\sum(\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i}))g(x_{i})\right|^{2}\leqslant\left(\frac{1}{N}\sum(\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i}))^{2}\right)\left(\frac{1}{N}\sum
g(x_{i})^{2}\right)$

By the strong law of large numbers, the second term $\frac{1}{N}\sum
g(x_{i})^{2}$ converges to $E_{P_{\theta}}g^{2}$ almost surely. So we have to
prove that the first term
$\frac{1}{N}\sum(\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i}))^{2}$ converges
to $0$ almost surely.

Since $w$ is bounded by assumption, we can write

$\displaystyle(\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i}))^{2}$
$\displaystyle\leqslant 2B\left|\mskip
1.0mu\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i})\right|$
$\displaystyle=2B\left|\mskip
1.0mu\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i})\right|_{+}+2B\left|\mskip
1.0mu\widehat{W}^{f}(x_{i})-W_{\theta}^{f}(x_{i})\right|_{-}$

where $B$ is the bound on $\left|\mskip 1.0muw\right|$. We will bound each of
these terms.

Let us abbreviate $q_{i}^{-}=\Pr_{x^{\prime}\sim
P_{\theta}}(f(x^{\prime})<f(x_{i}))$, $q_{i}^{+}=\Pr_{x^{\prime}\sim
P_{\theta}}(f(x^{\prime})\leqslant f(x_{i}))$,
$r_{i}^{-}=\\#\\{j\\!\leqslant\\!N,\,f(x_{j})<f(x_{i})\\}$,
$r_{i}^{+}=\\#\\{j\\!\leqslant\\!N,\,f(x_{j})\leqslant f(x_{i})\\}$.

By definition of $\widehat{W}^{f}$ we have

$\widehat{W}^{f}(x_{i})=\frac{1}{r_{i}^{+}-r_{i}^{-}}\sum_{k=r_{i}^{-}}^{r_{i}^{+}-1}w((k+1/2)/N)$

and moreover $W_{\theta}^{f}(x_{i})=w(q_{i}^{-})$ if $q_{i}^{-}=q_{i}^{+}$ or
$W_{\theta}^{f}(x_{i})=\frac{1}{q_{i}^{+}-q_{i}^{-}}\int_{q_{i}^{-}}^{q_{i}^{+}}w$
otherwise.

The Glivenko–Cantelli theorem states that $\sup_{i}\left|\mskip
1.0muq_{i}^{+}-r_{i}^{+}/N\right|$ tends to $0$ almost surely, and likewise
for $\sup_{i}\left|\mskip 1.0muq_{i}^{-}-r_{i}^{-}/N\right|$. So let $N$ be
large enough so that these errors are bounded by $\varepsilon$.

Since $w$ is non-increasing, we have $w(q_{i}^{-})\leqslant
w(r_{i}^{-}/N-\varepsilon)$. In the case $q_{i}^{-}\neq q_{i}^{+}$, we
decompose the interval $[q_{i}^{-};q_{i}^{+}]$ into $(r_{i}^{+}-r_{i}^{-})$
subintervals. The average of $w$ over each such subinterval is compared to a
term in the sum defining $w^{N}(x_{i})$: since $w$ is non-increasing, the
average of $w$ over the $k^{\text{th}}$ subinterval is at most
$w((r_{i}^{-}+k)/N-\varepsilon)$. So we get

$W_{\theta}^{f}(x_{i})\leqslant\frac{1}{r_{i}^{+}-r_{i}^{-}}\sum_{k=r_{i}^{-}}^{r_{i}^{+}-1}w(k/N-\varepsilon)$

so that

$W_{\theta}^{f}(x_{i})-\widehat{W}^{f}(x_{i})\leqslant\frac{1}{r_{i}^{+}-r_{i}^{-}}\sum_{k=r_{i}^{-}}^{r_{i}^{+}-1}(w(k/N-\varepsilon)-w((k+1/2)/N)).$

Let us sum over $i$, remembering that there are $(r_{i}^{+}-r_{i}^{-})$ values
of $j$ for which $f(x_{j})=f(x_{i})$. Taking the positive part, we get

$\frac{1}{N}\sum_{i=1}^{N}\left|\mskip
1.0muW_{\theta}^{f}(x_{i})-\widehat{W}^{f}(x_{i})\right|_{+}\leqslant\frac{1}{N}\sum_{k=0}^{N-1}(w(k/N-\varepsilon)-w((k+1/2)/N)).$

Since $w$ is non-increasing we have

$\frac{1}{N}\sum_{k=0}^{N-1}w(k/N-\varepsilon)\leqslant\int_{-\varepsilon-1/N}^{1-\varepsilon-1/N}w$

and

$\frac{1}{N}\sum_{k=0}^{N-1}w((k+1/2)/N)\geqslant\int_{1/2N}^{1+1/2N}w$

(we implicitly extend the range of $w$ so that $w(q)=w(0)$ for $q<0$). So we
have

$\frac{1}{N}\sum_{i=1}^{N}\left|\mskip
1.0muW_{\theta}^{f}(x_{i})-\widehat{W}^{f}(x_{i})\right|_{+}\leqslant\int_{-\varepsilon-1/N}^{1/2N}w-\int_{1-\varepsilon-1/N}^{1+1/2N}w\leqslant(2\varepsilon+3/N)B$

where $B$ is the bound on $\left|\mskip 1.0muw\right|$.

Reasoning symmetrically with $w(k/N+\varepsilon)$ and the inequalities
reversed, we get a similar bound for $\frac{1}{N}\sum\left|\mskip
1.0muW_{\theta}^{f}(x_{i})-\widehat{W}^{f}(x_{i})\right|_{-}$. This ends the
proof. ∎

### Proof of Proposition 5 (Quantile improvement)

Let us use the weight $w(u)=\mathbbm{1}_{u\leqslant q}$. Let $m$ be the value
of the $q$-quantile of $f$ under $P_{\theta^{t}}$. We want to show that the
value of the $q$-quantile of $f$ under $P_{\theta^{t+\hskip 0.25002pt\delta
t\hskip 0.25002pt}}$ is less than $m$, unless the gradient vanishes and the
IGO flow is stationary.

Let $p_{-}=\Pr_{x\sim P_{\theta^{t}}}(f(x)<m)$, $p_{m}=\Pr_{x\sim
P_{\theta^{t}}}(f(x)=m)$ and $p_{+}=\Pr_{x\sim P_{\theta^{t}}}(f(x)>m)$. By
definition of the quantile value we have $p_{-}+p_{m}\geqslant q$ and
$p_{+}+p_{m}\geqslant 1-q$. Let us assume that we are in the more complicated
case $p_{m}\neq 0$ (for the case $p_{m}=0$, simply remove the corresponding
terms).

We have $W_{\theta_{t}}^{f}(x)=1$ if $f(x)<m$, $W_{\theta_{t}}^{f}(x)=0$ if
$f(x)>m$ and
$W_{\theta_{t}}^{f}(x)=\frac{1}{p_{m}}\int_{p_{-}}^{p_{-}+p_{m}}w(u){\mathrm{d}}u=\frac{q-p_{-}}{p_{m}}$
if $f(x)=m$.

Using the same notation as above, let $g_{t}(\theta)=\int
W_{\theta^{t}}^{f}(x)\,P_{\theta}({\mathrm{d}}x)$. Decomposing this integral
on the three sets $f(x)<m$, $f(x)=m$ and $f(x)>m$, we get that
$g_{t}(\theta)=\Pr_{x\sim P_{\theta}}(f(x)<m)+\Pr_{x\sim
P_{\theta}}(f(x)=m)\frac{q-p_{-}}{p_{m}}$. In particular,
$g_{t}(\theta^{t})=q$.

Since we follow a gradient ascent of $g_{t}$, for $\hskip 0.50003pt\delta
t\hskip 0.50003pt$ small enough we have $g_{t}(\theta^{t+\hskip
0.35002pt\delta t\hskip 0.35002pt})>g_{t}(\theta^{t})$ unless the gradient
vanishes. If the gradient vanishes we have $\theta^{t+\hskip 0.35002pt\delta
t\hskip 0.35002pt}=\theta^{t}$ and the quantiles are the same. Otherwise we
get $g_{t}(\theta^{t+\hskip 0.35002pt\delta t\hskip
0.35002pt})>g_{t}(\theta^{t})=q$.

Since $\frac{q-p_{-}}{p_{m}}\leqslant\frac{(p_{-}+p_{m})-p_{-}}{p_{m}}=1$, we
have $g_{t}(\theta)\leqslant\Pr_{x\sim P_{\theta}}(f(x)<m)+\Pr_{x\sim
P_{\theta}}(f(x)=m)=\Pr_{x\sim P_{\theta}}(f(x)\leqslant m)$.

So $\Pr_{x\sim P_{\theta^{t+\hskip 0.25002pt\delta t\hskip
0.25002pt}}}(f(x)\leqslant m)\geqslant g_{t}(\theta^{t+\hskip 0.35002pt\delta
t\hskip 0.35002pt})>q$. This implies that, by definition, the $q$-quantile
value of $P_{\theta^{t+\hskip 0.25002pt\delta t\hskip 0.25002pt}}$ is smaller
than $m$.

### Proof of Proposition 10 (Speed of the IGO flow)

######  Lemma 19.

Let $X$ be a centered $L^{2}$ random variable with values in
${\mathbb{R}}^{d}$ and let $A$ be a real-valued $L^{2}$ random variable. Then

$\|\mathbb{E}(AX)\|\leqslant\sqrt{\lambda\operatorname{Var}A}$

where $\lambda$ is the largest eigenvalue of the covariance matrix of $X$
expressed in an orthonormal basis.

###### Proof of the lemma.

Let $v$ be any vector in ${\mathbb{R}}^{d}$; its norm satisfies

$\|v\|=\sup_{w,\,\|w\|\leqslant 1}(v\cdot w)$

and in particular

$\displaystyle\|\mathbb{E}(AX)\|$ $\displaystyle=\sup_{w,\,\|w\|\leqslant
1}(w\cdot\mathbb{E}(AX))$ $\displaystyle=\sup_{w,\,\|w\|\leqslant
1}\mathbb{E}(A\,(w\cdot X))$ $\displaystyle=\sup_{w,\,\|w\|\leqslant
1}\mathbb{E}((A-\mathbb{E}A)\,(w\cdot X))\quad\text{since $(w\cdot X)$ is
centered}$ $\displaystyle\leqslant\sup_{w,\,\|w\|\leqslant
1}\sqrt{\operatorname{Var}A}\,\sqrt{\mathbb{E}((w\cdot X)^{2})}$

by the Cauchy–Schwarz inequality and using the fact that $A$ is centered.

Now, in an orthonormal basis we have

$(w\cdot X)=\sum_{i}w_{i}X_{i}$

so that

$\displaystyle\mathbb{E}((w\cdot X)^{2})$
$\displaystyle=\mathbb{E}\left(({\textstyle\sum}_{i}w_{i}X_{i})({\textstyle\sum}_{j}w_{j}X_{j})\right)$
$\displaystyle={\textstyle\sum}_{i}{\textstyle\sum}_{j}\mathbb{E}(w_{i}X_{i}w_{j}X_{j})$
$\displaystyle={\textstyle\sum}_{i}{\textstyle\sum}_{j}w_{i}w_{j}\mathbb{E}(X_{i}X_{j})$
$\displaystyle={\textstyle\sum}_{i}{\textstyle\sum}_{j}w_{i}w_{j}C_{ij}$

with $C_{ij}$ the covariance matrix of $X$. The latter expression is the
scalar product $(w\cdot Cw)$. Since $C$ is a symmetric positive-semidefinite
matrix, $(w\cdot Cw)$ is at most $\lambda\|w\|^{2}$ with $\lambda$ the largest
eigenvalue of $C$. ∎

For the IGO flow we have
$\frac{{\mathrm{d}}\theta^{t}}{{\mathrm{d}}t}=\mathbb{E}_{x\sim
P_{\theta}}W_{\theta}^{f}(x)\tilde{\nabla}_{\theta}\ln P_{\theta}(x)$.

So applying the lemma, we get that the norm of
$\frac{{\mathrm{d}}\theta}{{\mathrm{d}}t}$ is at most
$\sqrt{\lambda\operatorname{Var}_{x\sim P_{\theta}}W_{\theta}^{f}(x)}$ where
$\lambda$ is the largest eivengalue of the covariance matrix of
$\tilde{\nabla}_{\theta}\ln P_{\theta}(x)$ (expressed in a coordinate system
where the Fisher matrix at the current point $\theta$ is the identity).

By construction of the quantiles, we have $\operatorname{Var}_{x\sim
P_{\theta}}W_{\theta}^{f}(x)\leqslant\operatorname{Var}_{[0,1]}w$ (with
equality unless there are ties). Indeed, for a given $x$, let $\mathcal{U}$ be
a uniform random variable in $[0,1]$ independent from $x$ and define the
random variable $Q=q^{-}(x)+(q^{+}(x)-q^{-}(x))\mathcal{U}$. Then $Q$ is
uniformly distributed between the upper and lower quantiles $q^{+}(x)$ and
$q^{-}(x)$ and thus we can rewrite $W_{\theta}^{f}(x)$ as
$\mathbb{E}(w(Q)|x)$. By the Jensen inequality we have
$\operatorname{Var}W_{\theta}^{f}(x)=\operatorname{Var}\mathbb{E}(w(Q)|x)\leqslant\operatorname{Var}w(Q)$.
In addition when $x$ is taken under $P_{\theta}$, $Q$ is uniformly distributed
in $[0,1]$ and thus $\operatorname{Var}w(Q)=\operatorname{Var}_{[0,1]}w$, i.e.
$\operatorname{Var}_{x\sim
P_{\theta}}W_{\theta}^{f}(x)\leqslant\operatorname{Var}_{[0,1]}w$.

Besides, consider the tangent space in $\Theta$-space at point $\theta^{t}$,
and let us choose an orthonormal basis in this tangent space for the Fisher
metric. Then, in this basis we have $\tilde{\nabla}_{i}\ln
P_{\theta}(x)=\partial_{i}\ln P_{\theta}(x)$. So the covariance matrix of
$\tilde{\nabla}\ln P_{\theta}(x)$ is $\mathbb{E}_{x\sim
P_{\theta}}(\partial_{i}\ln P_{\theta}(x)\partial_{j}\ln P_{\theta}(x))$,
which is equal to the Fisher matrix by definition. So this covariance matrix
is the identity, whose largest eigenvalue is $1$.

### Proof of Proposition 12 (Noisy IGO)

On the one hand, let $P_{\theta}$ be a family of distributions on $X$. The IGO
algorithm (13) applied to a random function $f(x)=\tilde{f}(x,\omega)$ where
$\omega$ is a random variable uniformly distributed in $[0,1]$ reads

$\theta^{t+\hskip 0.35002pt\delta t\hskip 0.35002pt}=\theta^{t}+\hskip
0.50003pt\delta t\hskip
0.50003pt\sum_{i=1}^{N}\hat{w_{i}}\widetilde{\nabla}_{\theta}\ln
P_{\theta}(x_{i})$ (60)

where $x_{i}\sim P_{\theta}$ and $\hat{w_{i}}$ is according to (12) where
ranking is applied to the values $\tilde{f}(x_{i},\omega_{i})$, with
$\omega_{i}$ uniform variables in $[0,1]$ independent from $x_{i}$ and from
each other.

On the other hand, for the IGO algorithm using $P_{\theta}\otimes U_{[0,1]}$
and applied to the deterministic function $\tilde{f}$, $\hat{w_{i}}$ is
computed using the ranking according to the $\tilde{f}$ values of the sampled
points $\tilde{x}_{i}=(x_{i},\omega_{i})$, and thus coincides with the one in
(60).

Besides,

$\widetilde{\nabla}_{\theta}\ln P_{\theta\otimes
U_{[0,1]}}(\tilde{x}_{i})=\widetilde{\nabla}_{\theta}\ln
P_{\theta}(x_{i})+\underbrace{\widetilde{\nabla}_{\theta}\ln
U_{[0,1]}(\omega_{i})}_{=0}$

and thus the IGO algorithm update on space $X\times[0,1]$, using the family of
distributions $\tilde{P}_{\theta}=P_{\theta}\otimes U_{[0,1]}$, applied to the
deterministic function $\tilde{f}$, coincides with (60).

### Proof of Theorem 13 (Natural gradient as ML with infinitesimal weights)

We begin with a calculus lemma (proof omitted).

######  Lemma 20.

Let $f$ be real-valued function on a finite-dimensional vector space $E$
equipped with a definite positive quadratic form $\|\cdot\|^{2}$. Assume $f$
is smooth and has at most quadratic growth at infinity. Then, for any $x\in
E$, we have

$\nabla f(x)=\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}\operatorname*{arg\,max}_{h}\left\\{f(x+h)-\frac{1}{2\varepsilon}\|h\|^{2}\right\\}$

where $\nabla$ is the gradient associated with the norm $\|\cdot\|$.
Equivalently,

$\operatorname*{arg\,max}_{y}\left\\{f(y)-\frac{1}{2\varepsilon}\|y-x\|^{2}\right\\}=x+\varepsilon\nabla
f(x)+O(\varepsilon^{2})$

when $\varepsilon\to 0$.

We are now ready to prove Theorem 13. Let $W$ be a function of $x$, and fix
some $\theta_{0}$ in $\Theta$.

We need some regularity assumptions: we assume that the parameter space
$\Theta$ is non-degenerate (no two points $\theta\in\Theta$ define the same
probability distribution) and proper (the map $P_{\theta}\mapsto\theta$ is
continuous). We also assume that the map $\theta\mapsto P_{\theta}$ is smooth
enough, so that $\int\log P_{\theta}(x)\,W(x)\,P_{\theta_{0}}({\mathrm{d}}x)$
is a smooth function of $\theta$. (These are restrictions on
$\theta$-regularity: this does not mean that $W$ has to be continuous as a
function of $x$.)

The two statements of Theorem 13 using a sum and an integral have similar
proofs, so we only include the first. For $\varepsilon>0$, let $\theta$ be the
solution of

$\theta=\operatorname*{arg\,max}\Bigg{\\{}(1-\varepsilon)\int\log
P_{\theta}(x)\,P_{\theta_{0}}({\mathrm{d}}x)+\varepsilon\int\log
P_{\theta}(x)\,W(x)\,P_{\theta_{0}}({\mathrm{d}}x)\Bigg{\\}}.$

Then we have

$\displaystyle\theta$
$\displaystyle=\operatorname*{arg\,max}\Bigg{\\{}\int\log
P_{\theta}(x)\,P_{\theta_{0}}({\mathrm{d}}x)+\varepsilon\int\log
P_{\theta}(x)\,(W(x)-1)\,P_{\theta_{0}}({\mathrm{d}}x)\Bigg{\\}}$
$\displaystyle=\operatorname*{arg\,max}\Bigg{\\{}\int\log
P_{\theta}(x)\,P_{\theta_{0}}({\mathrm{d}}x)-\int\log
P_{\theta_{0}}(x)\,P_{\theta_{0}}({\mathrm{d}}x)+\varepsilon\int\log
P_{\theta}(x)\,(W(x)-1)\,P_{\theta_{0}}({\mathrm{d}}x)\Bigg{\\}}$ (because the
added term does not depend on $\theta$)
$\displaystyle=\operatorname*{arg\,max}\Bigg{\\{}-\mathrm{KL}\\!\left(P_{\theta_{0}}\,||\,P_{\theta}\right)+\varepsilon\int\log
P_{\theta}(x)\,(W(x)-1)\,P_{\theta_{0}}({\mathrm{d}}x)\Bigg{\\}}\ $
$\displaystyle=\operatorname*{arg\,max}\Bigg{\\{}-\frac{1}{\varepsilon}\mathrm{KL}\\!\left(P_{\theta_{0}}\,||\,P_{\theta}\right)+\int\log
P_{\theta}(x)\,(W(x)-1)\,P_{\theta_{0}}({\mathrm{d}}x)\Bigg{\\}}$

When $\varepsilon\to 0$, the first term exceeds the second one if
$\mathrm{KL}\\!\left(P_{\theta_{0}}\,||\,P_{\theta}\right)$ is too large
(because $W$ is bounded), and so
$\mathrm{KL}\\!\left(P_{\theta_{0}}\,||\,P_{\theta}\right)$ tends to $0$. So
we can assume that $\theta$ is close to $\theta_{0}$.

When $\theta=\theta_{0}+\delta\theta$ is close to $\theta_{0}$, we have

$\mathrm{KL}\\!\left(P_{\theta_{0}}\,||\,P_{\theta}\right)=\frac{1}{2}\sum
I_{ij}(\theta_{0})\,\delta\theta_{i}\,\delta\theta_{j}+O(\delta\theta^{3})$

with $I_{ij}(\theta_{0})$ the Fisher matrix at $\theta_{0}$. (This actually
holds both for $\mathrm{KL}\\!\left(P_{\theta_{0}}\,||\,P_{\theta}\right)$ and
$\mathrm{KL}\\!\left(P_{\theta}\,||\,P_{\theta_{0}}\right)$.)

Thus, we can apply the lemma above using the Fisher metric $\sum
I_{ij}(\theta_{0})\,\delta\theta_{i}\,\delta\theta_{j}$, and working on a
small neighborhood of $\theta_{0}$ in $\theta$-space (which can be identified
with ${\mathbb{R}}^{\dim\Theta}$). The lemma states that the argmax above is
attained at

$\theta=\theta_{0}+\varepsilon\widetilde{\nabla}_{\theta}\int\log
P_{\theta}(x)\,(W(x)-1)\,P_{\theta_{0}}({\mathrm{d}}x)$

up to $O(\varepsilon^{2})$, with $\widetilde{\nabla}$ the natural gradient.

Finally, the gradient cancels the constant $-1$ because
$\int(\widetilde{\nabla}\log P_{\theta})\,P_{\theta_{0}}=0$ at
$\theta=\theta_{0}$. This proves Theorem 13.

### Proof of Theorem 15 (IGO, CEM and IGO-ML)

Let $P_{\theta}$ be a family of probability distributions of the form

$P_{\theta}(x)=\frac{1}{Z(\theta)}\exp\left(\sum\theta_{i}T_{i}(x)\right)\,H({\mathrm{d}}x)$

where $T_{1},\ldots,T_{k}$ is a finite family of functions on $X$ and
$H({\mathrm{d}}x)$ is some reference measure on $X$. We assume that the family
of functions $(T_{i})_{i}$ together with the constant function $T_{0}(x)=1$,
are linearly independent. This prevents redundant parametrizations where two
values of $\theta$ describe the same distribution; this also ensures that the
Fisher matrix $\operatorname{Cov}(T_{i},T_{j})$ is invertible.

The IGO update (14) in the parametrization $\bar{T}_{i}$ is a sum of terms of
the form

$\widetilde{\nabla}_{\bar{T}_{i}}\ln P(x).$

So we will compute the natural gradient $\widetilde{\nabla}_{\bar{T}_{i}}$ in
those coordinates. We first need some general results about the Fisher metric
for exponential families.

The next proposition gives the expression of the Fisher scalar product between
two tangent vectors $\delta P$ and $\delta^{\prime}P$ of a statistical
manifold of exponential distributions. It is one way to express the duality
between the coordinates $\theta_{i}$ and $\bar{T}_{i}$ (compare [AN00, (3.30)
and Section 3.5]).

######  Proposition 21.

Let $\delta\theta_{i}$ and $\delta^{\prime}\theta_{i}$ be two small variations
of the parameters $\theta_{i}$. Let $\delta P(x)$ and $\delta^{\prime}P(x)$ be
the resulting variations of the probability distribution $P$, and
$\delta\bar{T}_{i}$ and $\delta^{\prime}\bar{T}_{i}$ the resulting variations
of $\bar{T}_{i}$. Then the scalar product, in Fisher information metric,
between the tangent vectors $\delta P$ and $\delta^{\prime}P$, is

$\langle\delta
P,\delta^{\prime}P\rangle=\sum_{i}\delta\theta_{i}\,\delta^{\prime}\bar{T}_{i}=\sum_{i}\delta^{\prime}\theta_{i}\,\delta\bar{T}_{i}.$

###### Proof.

By definition of the Fisher metric:

$\displaystyle\langle\delta P,\delta^{\prime}P\rangle$
$\displaystyle=\sum_{i,j}I_{ij}\,\delta\theta_{i}\,\delta^{\prime}\theta_{j}$
$\displaystyle=\sum_{i,j}\delta\theta_{i}\,\delta^{\prime}\theta_{j}\int_{x}\frac{\partial\ln
P(x)}{\partial\theta_{i}}\,\frac{\partial\ln P(x)}{\partial\theta_{j}}\,P(x)$
$\displaystyle=\int_{x}\sum_{i}\frac{\partial\ln
P(x)}{\partial\theta_{i}}\delta\theta_{i}\,\sum_{j}\frac{\partial\ln
P(x)}{\partial\theta_{j}}\delta^{\prime}\theta_{j}\,P(x)$
$\displaystyle=\int_{x}\sum_{i}\frac{\partial\ln
P(x)}{\partial\theta_{i}}\,\delta\theta_{i}\;\delta^{\prime}(\ln P(x))\,P(x)$
$\displaystyle=\int_{x}\sum_{i}\frac{\partial\ln
P(x)}{\partial\theta_{i}}\,\delta\theta_{i}\;\delta^{\prime}P(x)$
$\displaystyle=\int_{x}\sum_{i}(T_{i}(x)-\bar{T}_{i})\delta\theta_{i}\;\delta^{\prime}P(x)\qquad\text{by
\eqref{eq:gradexp}}$
$\displaystyle=\sum_{i}\delta\theta_{i}\left(\int_{x}T_{i}(x)\,\delta^{\prime}P(x)\right)-\sum_{i}\delta\theta_{i}\bar{T}_{i}\int_{x}\delta^{\prime}P(x)$
$\displaystyle=\sum_{i}\delta\theta_{i}\,\delta^{\prime}\bar{T}_{i}$

because $\int_{x}\delta^{\prime}P(x)=0$ since the total mass of $P$ is $1$,
and $\int_{x}T_{i}(x)\,\delta^{\prime}P(x)=\delta^{\prime}\bar{T}_{i}$ by
definition of $\bar{T}_{i}$. ∎

######  Proposition 22.

Let $f$ be a function on the statistical manifold of an exponential family as
above. Then the components of the natural gradient w.r.t. the expectation
parameters are given by the vanilla gradient w.r.t. the natural parameters:

$\widetilde{\nabla}_{\bar{T}_{i}}f=\frac{\partial f}{\partial\theta_{i}}$

and conversely

$\widetilde{\nabla}_{\theta_{i}}f=\frac{\partial f}{\partial\bar{T}_{i}}.$

(Beware this does _not_ mean that the gradient ascent in any of those
parametrizations is the vanilla gradient ascent.)

We could not find a reference for this result, though we think it known.

###### Proof.

By definition, the natural gradient $\widetilde{\nabla}f$ of a function $f$ is
the unique tangent vector $\delta P$ such that that, for any other tangent
vector $\delta^{\prime}P$, we have

$\delta^{\prime}f=\langle\delta P,\delta^{\prime}P\rangle$

with $\langle\cdot,\cdot\rangle$ the scalar product associated with the Fisher
metric. We want to compute this natural gradient in coordinates $\bar{T}_{i}$,
so we are interested in the variations $\delta\bar{T}_{i}$ associated with
$\delta P$.

By Proposition 21, the scalar product above is

$\langle\delta
P,\delta^{\prime}P\rangle=\sum\delta\bar{T}_{i}\,\delta^{\prime}\theta_{i}$

where $\delta\bar{T}_{i}$ is the variation of $\bar{T}_{i}$ associated with
$\delta P$, and $\delta^{\prime}\theta_{i}$ the variation of $\theta_{i}$
associated with $\delta^{\prime}P$.

On the other hand we have $\delta^{\prime}f=\sum_{i}\frac{\partial
f}{\partial\theta_{i}}\,\delta^{\prime}\theta_{i}$. So we must have

$\sum_{i}\delta\bar{T}_{i}\,\delta^{\prime}\theta_{i}=\sum_{i}\frac{\partial
f}{\partial\theta_{i}}\,\delta^{\prime}\theta_{i}$

for any $\delta^{\prime}P$, which leads to

$\delta\bar{T}_{i}=\frac{\partial f}{\partial\theta_{i}}$

as needed. The converse relation is proved _mutatis mutandis_. ∎

Back to the proof of Theorem 15. We can now compute the desired terms:

$\widetilde{\nabla}_{\bar{T}_{i}}\ln P(x)=\frac{\partial\ln
P(x)}{\partial\theta_{i}}=T_{i}(x)-\bar{T}_{i}$

by (16). This proves the first statement (27) in Theorem 15 about the form of
the IGO update in these parameters.

The other statements follow easily from this together with the additional fact
(26) that, for any set of weights $a_{i}$ with $\sum a_{i}=1$, the value
$T^{\ast}=\sum_{i}a(i)T(x_{i})$ is the maximum likelihood estimate of
$\sum_{i}a(i)\log P(x_{i})$.

## References

  * [AHS85] D.H. Ackley, G.E. Hinton, and T.J. Sejnowski, _A learning algorithm for Boltzmann machines_ , Cognitive Science 9 (1985), no. 1, 147–169.
  * [Ama98] Shun-Ichi Amari, _Natural gradient works efficiently in learning_ , Neural Comput. 10 (1998), 251–276.
  * [AN00] Shun-ichi Amari and Hiroshi Nagaoka, _Methods of information geometry_ , Translations of Mathematical Monographs, vol. 191, American Mathematical Society, Providence, RI, 2000, Translated from the 1993 Japanese original by Daishi Harada.
  * [ANOK10] Youhei Akimoto, Yuichi Nagata, Isao Ono, and Shigenobu Kobayashi, _Bidirectional relation between CMA evolution strategies and natural evolution strategies_ , Proceedings of Parallel Problem Solving from Nature - PPSN XI, Lecture Notes in Computer Science, vol. 6238, Springer, 2010, pp. 154–163.
  * [AO08] Ravi P. Agarwal and Donal O’Regan, _An introduction to ordinary differential equations_ , Springer, 2008.
  * [Arn06] D.V. Arnold, _Weighted multirecombination evolution strategies_ , Theoretical computer science 361 (2006), no. 1, 18–37.
  * [Bal94] Shumeet Baluja, _Population based incremental learning: A method for integrating genetic search based function optimization and competitve learning_ , Tech. Report CMU-CS-94-163, Carnegie Mellon Report, 1994.
  * [BC95] Shumeet Baluja and Rich Caruana, _Removing the genetics from the standard genetic algorithm_ , Proceedings of ICML’95, 1995, pp. 38–46.
  * [Ber00] A. Berny, _Selection and reinforcement learning for combinatorial optimization_ , Parallel Problem Solving from Nature PPSN VI (Marc Schoenauer, Kalyanmoy Deb, Günther Rudolph, Xin Yao, Evelyne Lutton, Juan Merelo, and Hans-Paul Schwefel, eds.), Lecture Notes in Computer Science, vol. 1917, Springer Berlin Heidelberg, 2000, pp. 601–610.
  * [Ber02] Arnaud Berny, _Boltzmann machine for population-based incremental learning_ , ECAI, 2002, pp. 198–202.
  * [Bey01] H.-G. Beyer, _The theory of evolution strategies_ , Natural Computing Series, Springer-Verlag, 2001.
  * [BLS07] J. Branke, C. Lode, and J.L. Shapiro, _Addressing sampling errors and diversity loss in umda_ , Proceedings of the 9th annual conference on Genetic and evolutionary computation, ACM, 2007, pp. 508–515.
  * [BS02] H.G. Beyer and H.P. Schwefel, _Evolution strategies—a comprehensive introduction_ , Natural computing 1 (2002), no. 1, 3–52.
  * [CT06] Thomas M. Cover and Joy A. Thomas, _Elements of information theory_ , second ed., Wiley-Interscience [John Wiley & Sons], Hoboken, NJ, 2006.
  * [dBKMR05] Pieter-Tjerk de Boer, Dirk P. Kroese, Shie Mannor, and Reuven Y. Rubinstein, _A tutorial on the cross-entropy method_ , Annals OR 134 (2005), no. 1, 19–67.
  * [Gha04] Zoubin Ghahramani, _Unsupervised learning_ , Advanced Lectures on Machine Learning (Olivier Bousquet, Ulrike von Luxburg, and Gunnar Rätsch, eds.), Lecture Notes in Computer Science, vol. 3176, Springer Berlin / Heidelberg, 2004, pp. 72–112.
  * [GSS+10] Tobias Glasmachers, Tom Schaul, Yi Sun, Daan Wierstra, and Jürgen Schmidhuber, _Exponential natural evolution strategies_ , GECCO, 2010, pp. 393–400.
  * [Han06a] N. Hansen, _An analysis of mutative $\sigma$-self-adaptation on linear fitness functions_, Evolutionary Computation 14 (2006), no. 3, 255–275.
  * [Han06b] , _The CMA evolution strategy: a comparing review_ , Towards a new evolutionary computation. Advances on estimation of distribution algorithms (J.A. Lozano, P. Larranaga, I. Inza, and E. Bengoetxea, eds.), Springer, 2006, pp. 75–102.
  * [Han09] , _Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed_ , Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers (New York, NY, USA), GECCO ’09, ACM, 2009, pp. 2389–2396.
  * [Hin02] G.E. Hinton., _Training products of experts by minimizing contrastive divergence_ , Neural Computation 14 (2002), 1771–1800.
  * [HJ61] R. Hooke and T.A. Jeeves, _“direct search” solution of numerical and statistical problems_ , Journal of the ACM 8 (1961), 212–229.
  * [HK04] N. Hansen and S. Kern, _Evaluating the CMA evolution strategy on multimodal test functions_ , Parallel Problem Solving from Nature PPSN VIII (X. Yao et al., eds.), LNCS, vol. 3242, Springer, 2004, pp. 282–291.
  * [HMK03] N. Hansen, S.D. Müller, and P. Koumoutsakos, _Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)_ , Evolutionary Computation 11 (2003), no. 1, 1–18.
  * [HO96] N. Hansen and A. Ostermeier, _Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation_ , ICEC96, IEEE Press, 1996, pp. 312–317.
  * [HO01] Nikolaus Hansen and Andreas Ostermeier, _Completely derandomized self-adaptation in evolution strategies_ , Evolutionary Computation 9 (2001), no. 2, 159–195.
  * [JA06] G.A. Jastrebski and D.V. Arnold, _Improving evolution strategies through active covariance matrix adaptation_ , Evolutionary Computation, 2006. CEC 2006\. IEEE Congress on, IEEE, 2006, pp. 2814–2821.
  * [Jef46] Harold Jeffreys, _An invariant form for the prior probability in estimation problems_ , Proc. Roy. Soc. London. Ser. A. 186 (1946), 453–461.
  * [Kul97] Solomon Kullback, _Information theory and statistics_ , Dover Publications Inc., Mineola, NY, 1997, Reprint of the second (1968) edition.
  * [LL02] P. Larranaga and J.A. Lozano, _Estimation of distribution algorithms: A new tool for evolutionary computation_ , Springer Netherlands, 2002.
  * [LRMB07] Nicolas Le Roux, Pierre-Antoine Manzagol, and Yoshua Bengio, _Topmoumoute online natural gradient algorithm_ , NIPS, 2007.
  * [MMS08] Luigi Malagò, Matteo Matteucci, and Bernardo Dal Seno, _An information geometry perspective on estimation of distribution algorithms: boundary analysis_ , GECCO (Companion), 2008, pp. 2081–2088.
  * [NM65] John Ashworth Nelder and R Mead, _A simplex method for function minimization_ , The Computer Journal (1965), 308–313.
  * [PGL02] M. Pelikan, D.E. Goldberg, and F.G. Lobo, _A survey of optimization by building and using probabilistic models_ , Computational optimization and applications 21 (2002), no. 1, 5–20.
  * [Rao45] Calyampudi Radhakrishna Rao, _Information and the accuracy attainable in the estimation of statistical parameters_ , Bull. Calcutta Math. Soc. 37 (1945), 81–91.
  * [Rec73] I. Rechenberg, _Evolutionstrategie: Optimierung technisher systeme nach prinzipien des biologischen evolution_ , Frommann-Holzboog Verlag, Stuttgart, 1973\.
  * [RK04] R.Y. Rubinstein and D.P. Kroese, _The cross-entropy method: a unified approach to combinatorial optimization, monte-carlo simulation, and machine learning_ , Springer-Verlag New York Inc, 2004.
  * [Rub99] Reuven Rubinstein, _The cross-entropy method for combinatorial and continuous optimization_ , Methodology and Computing in Applied Probability 1 (1999), 127–190, 10.1023/A:1010091220143.
  * [SA90] F. Silva and L. Almeida, _Acceleration techniques for the backpropagation algorithm_ , Neural Networks (1990), 110–119.
  * [Sch92] Laurent Schwartz, _Analyse. II_ , Collection Enseignement des Sciences [Collection: The Teaching of Science], vol. 43, Hermann, Paris, 1992, Calcul différentiel et équations différentielles, With the collaboration of K. Zizi.
  * [Sch95] H.-P. Schwefel, _Evolution and optimum seeking_ , Sixth-generation computer technology series, John Wiley & Sons, Inc. New York, NY, USA, 1995.
  * [SHI09] Thorsten Suttorp, Nikolaus Hansen, and Christian Igel, _Efficient covariance matrix update for variable metric evolution strategies_ , Machine Learning 75 (2009), no. 2, 167–197.
  * [Smo86] P. Smolensky, _Information processing in dynamical systems: foundations of harmony theory_ , Parallel Distributed Processing (D. Rumelhart and J. McClelland, eds.), vol. 1, MIT Press, Cambridge, MA, USA, 1986, pp. 194–281.
  * [SWSS09] Yi Sun, Daan Wierstra, Tom Schaul, and Juergen Schmidhuber, _Efficient natural evolution strategies_ , Proceedings of the 11th Annual conference on Genetic and evolutionary computation (New York, NY, USA), GECCO ’09, ACM, 2009, pp. 539–546.
  * [Tho00] H. Thorisson, _Coupling, stationarity, and regeneration_ , Springer, 2000.
  * [Tor97] Virginia Torczon, _On the convergence of pattern search algorithms_ , SIAM Journal on optimization 7 (1997), no. 1, 1–25.
  * [WAS04] Michael Wagner, Anne Auger, and Marc Schoenauer, _EEDA : A New Robust Estimation of Distribution Algorithms_ , Research Report RR-5190, INRIA, 2004\.
  * [Whi89] D. Whitley, _The genitor algorithm and selection pressure: Why rank-based allocation of reproductive trials is best_ , Proceedings of the third international conference on Genetic algorithms, 1989, pp. 116–121.
  * [WSPS08] Daan Wierstra, Tom Schaul, Jan Peters, and Jürgen Schmidhuber, _Natural evolution strategies_ , IEEE Congress on Evolutionary Computation, 2008, pp. 3381–3387.