Black-Box Optimization Using Geodesics in Statistical Manifolds ^†

Abstract

Information geometric optimization (IGO) is a general framework for stochastic optimization problems aiming at limiting the influence of arbitrary parametrization choices: the initial problem is transformed into the optimization of a smooth function on a Riemannian manifold, defining a parametrization-invariant first order differential equation and, thus, yielding an approximately parametrization-invariant algorithm (up to second order in the step size). We define the geodesic IGO update, a fully parametrization-invariant algorithm using the Riemannian structure, and we compute it for the manifold of Gaussians, thanks to Noether’s theorem. However, in similar algorithms, such as CMA-ES (Covariance Matrix Adaptation - Evolution Strategy) and xNES (exponential Natural Evolution Strategy), the time steps for the mean and the covariance are decoupled. We suggest two ways of doing so: twisted geodesic IGO (GIGO) and blockwise GIGO. Finally, we show that while the xNES algorithm is not GIGO, it is an instance of blockwise GIGO applied to the mean and covariance matrix separately. Therefore, xNES has an almost parametrization-invariant description.

1. Introduction

Consider an objective function f: X → ℝ to be minimized. We suppose we have absolutely no knowledge about f: the only thing we can do is ask for its value at any point x ∈ X (black-box optimization) and that the evaluation of f is a costly operation. We are going to study algorithms that can be described in the IGO framework (see [1]).

We consider the following optimization procedure:

We choose (P_θ)_θ∈Θ a family of probability distributions (which will be given a Riemannian manifold structure, following [2]) on X and an initial probability distribution $P_{{\theta }^0}$. Now, we replace f by F: Θ → ℝ (for example $F(\theta )=E_{x~P_{\theta }}[f(x)]$), and we optimize F by gradient descent, corresponding to the gradient flow:

$$\frac{\mathrm{d}{\theta }^t}{\mathrm{d}t}=-{\nabla }_{\theta }{E}_{x~P_{\theta }}[f(x)].$$

(1)

However, because of the gradient, this equation depends entirely on the parametrization we chose for Θ, which is disturbing: we do not want to have two different updates, because we chose different parameters to represent the objects with which we are working. Moreover, in the case of a function with several local minima, changing the parametrization can change the attained optimum (see [3], for example). That is why invariance is a design principle behind IGO. More precisely, we want invariance with respect to monotone transformations of f and invariance under reparametrization of Θ.

The IGO framework uses the geometry of the family Θ, which is given by the Fisher metric to provide a differential equation on θ with the desired properties, but because of the discretization of time needed to obtain an explicit algorithm, we lose invariance under reparametrization of θ: two IGO algorithms applied to the same function to be optimized, but with different parametrizations, coincide only at first order in the step size. A possible solution to this problem is geodesic IGO (GIGO), introduced here (see also IGO-Maximum Likelihoodin [1], for example.): the initial direction of the update at each step of the algorithm remains the same as in IGO, but instead of moving straight for the chosen parametrization, we use the Riemannian manifold structure of our family of probability distributions (see [2]) by following its geodesics.

Finding the geodesics of a Riemannian manifold is not always easy, but Noether’s theorem will allow us to obtain quantities that are preserved along the geodesics, thus allowing, in the case of Gaussian distributions, one to obtain a first order differential equation satisfied by the geodesics, which makes their computation easier.

Although the geodesic IGO algorithm is not, strictly speaking, parametrization-invariant when no closed form for the geodesics is known, it is possible to compute them at arbitrary precision without increasing the numbers of objective function calls.

The first two sections are preliminaries: in Section 2, we recall the IGO algorithm, introduced in [1], and in Section 3, after a reminder about Riemannian geometry, we state Noether’s theorem, which will be our main tool to compute the GIGO update for Gaussian distributions.

In Section 4, we consider Gaussian distributions with a covariance matrix proportional to the identity matrix: this space is isometric to the hyperbolic space, and the geodesics of the latter are known.

In Section 5.1, we consider the general Gaussian case, and we use Noether’s theorem to obtain two different sets of equations to compute the GIGO update. The equations are known (see [4–6]), but the connection with Noether’s theorem has not been mentioned. We then give the explicit solution for these equations, from [5].

In Section 6, we recall quickly the xNES and CMA-ESupdates, and we introduce a slight modification of the IGO algorithm to incorporate the direction-dependent learning rates used in CMA-ESand xNES. We then compare these different algorithms and prove that xNES is not GIGO in general, and we finally introduce a new family of algorithms extending GIGO and recovering xNES from abstract principles.

Finally, Section 7 presents numerical experiments, which suggest that when using GIGO with Gaussian distributions, the step size must be chosen carefully.

2. Definitions: IGO, GIGO

In this section, we recall what the IGO framework is and we define the geodesic IGO update. Consider again Equation (1):

$$\frac{\mathrm{d}{\theta }^t}{\mathrm{d}t}=-{\nabla }_{\theta }{E}_{x~P_{\theta }}[f(x)].$$

As we saw in the Introduction:

• The gradient depends on the parametrization of our space of probability distributions (see Section 2.3 for an example).
• The equation is not invariant under monotone transformations of f. For example, the optimization for 10f moves ten times faster than the optimization for f.

In this section, we recall how IGO deals with this (see [1] for a better presentation).

2.1. Invariance under Reparametrization of θ: Fisher Metric

In order to achieve invariance under reparametrization of θ, it is possible to turn our family of probability distributions into a Riemannian manifold (this is the main topic of information geometry; see [2]), which allows us to use a canonical, parametrization-invariant gradient (called the natural gradient).

Definition 1. Let P, Q be two probability distributions on X. The Kullback–Leibler divergence of Q from P is defined by:

$$\mathrm{KL}(Q\Vert P)={\displaystyle {\int }_X\ln (\frac{Q(x)}{P(x)})dQ(x)}.$$

(2)

By definition, it does not depend on the parametrization. It is not symmetrical, but if for all x, the application θ ⟼ P_θ(x) is C², then a second-order expansion yields:

$$\mathrm{KL}(P_{\theta +\mathrm{d}\theta }\Vert P_{\theta })=\frac{1}{2}{\displaystyle \sum _{i,j}{I}_{ij}(\theta )\mathrm{d}{\theta }_i\mathrm{d}{\theta }_j+o(\mathrm{d}{\theta }^2)},$$

(3)

where:

$$I_{ij}(\theta )={\displaystyle {\int }_X\frac{\partial \ln P_{\theta }(x)}{\partial {\theta }_i}\frac{\partial \ln P_{\theta }(x)}{\partial {\theta }_j}\mathrm{d}{P}_{\theta }(x)=-{\displaystyle {\int }_X\frac{{\partial }^2\ln P_{\theta }(x)}{\partial {\theta }_i\partial {\theta }_j}\mathrm{d}{P}_{\theta }(x)}.}$$

(4)

This is enough to endow the family (P_θ)_θ∈Θ with a Riemannian manifold structure: a Riemannian manifold M is a differentiable manifold, which can be seen as pieces of ℝⁿ glued together, with a metric. The metric at x is a symmetric positive-definite quadratic form on the tangent space of M at x: it indicates how expensive it is to move in a given direction on the manifold. We will think of the updates of the algorithms that we will be studying as paths on M.

The matrix I(θ) is called the “Fisher information matrix”, and the metric it defines is called the “Fisher metric”.

Given a metric, it is possible to define a gradient attached to this metric; the key property of the gradient is that for any smooth function f:

$$f\left(x+h\right)=f\left(x\right)+{\displaystyle \sum _i{h}_i\frac{\partial f}{\partial x_i}+o\left(\Vert h\Vert \right)}=f\left(x\right)+\langle x\rangle +\langle h,\nabla f\left(x\right)\rangle +o\left(\Vert h\Vert \right),$$

(5)

where $\langle x,y\rangle =x^T{I}_y$ is the dot product in metric I. Therefore, in order to keep the property of Equation (5), we must have $\nabla f=I^{-1}\frac{\partial f}{\partial x}$.

We have therefore the following gradient (called the “natural gradient”; see [2]):

$${\tilde{\nabla }}_{\theta }=I^{-1}(\theta )\frac{\partial }{\partial \theta },$$

(6)

and since the Kullback–Leibler divergence does not depend on the parametrization, neither does the natural gradient.

Later in this paper, we will study families of Gaussian distributions. The following proposition gives the Fisher metric for these families.

Proposition 1. Let (P_θ)_θ∈Θ be a family of normal probability distributions:$P_{\theta }=N(\mu (\theta ),\mathrm{\Sigma }(\theta )).$ If μ and Σ are C¹, the Fisher metric is given by:

$$I_{i,j}(\theta )=\frac{\partial {\mu }^T}{\partial {\theta }_i}{\mathrm{\Sigma }}^{-1}\frac{\partial \mu }{\partial {\theta }_j}+\frac{1}{2}\mathrm{tr}\left({\mathrm{\Sigma }}^{-1}\frac{\partial \mathrm{\Sigma }}{\partial {\theta }_i}{\mathrm{\Sigma }}^{-1}\frac{\partial \mathrm{\Sigma }}{\partial {\theta }_j}\right).$$

(7)

Proof. This is a non-trivial calculation. See [7] or [8] for more details.

As we will often be working with Gaussian distributions, we introduce the following notation:

Notation 1. ${\mathbb{G}}_d$ is the manifold of Gaussian distributions in dimension d, equipped with the Fisher metric.${\tilde{\mathbb{G}}}_d$ is the manifold of Gaussian distributions in dimension d, with the covariance matrix proportional to identity in the canonical basis of ℝ^d, equipped with the Fisher metric.

2.2. IGO Flow, IGO Algorithm

In IGO [1], invariance with respect to monotone transformations is achieved by replacing f by the following transform; we set:

$$q(x)=P_{x^{\prime }~P_{\theta }}(f(x^{\prime })\le f(x)),$$

(8)

a non-increasing function w: [0; 1] → ℝ is chosen (the selection scheme), and finally, $W_{\theta }^f(x)=w(q(x))$ (this definition has to be slightly changed if the probability of a tie is not zero, see [1] for more details). By performing a gradient descent on $E_{x~P_{\theta }}\left[W_{{\theta }^t}^f(x)\right]$, we obtain the “IGO flow”:

$$\frac{\mathrm{d}{\theta }^t}{\mathrm{d}t}={\tilde{\nabla }}_{\theta }{\displaystyle {\int }_X{W}_{{\theta }^t}^f(x)}{P}_{\theta }(\mathrm{d}x)={\displaystyle {\int }_X{W}_{{\theta }^t}^f(x)}{\tilde{\nabla }}_{\theta }\ln P_{\theta }(x)P_{{\theta }^t}(\mathrm{d}x).$$

(9)

Notice that the function we are optimizing is $E_{x~P_{\theta }}\left[W_{{\theta }^t}^f(x)\right]$ and not $E_{x~P_{\theta }}\left[W_{\theta }^f(x)\right]$ (the second function is constant and always equal to ${\int }_0^{1}w$). In particular, the function for which we are performing the gradient descent changes at each step, although their optimum (a Dirac at the minimum of f) does not: the IGO flow is not a gradient flow; it is simply a vector flow given by the gradient of interrelated functions.

For practical implementation, the integral in (9) has to be approximated. For the integral itself, the Monte-Carlo method is used; N values (x₁, …, x_N) are sampled from the distribution $P_{{\theta }^t}$, and the integral becomes:

$$\frac{1}{N}{\displaystyle \sum _{i=1}^N{W}_{{\theta }^t}^f(x_i){\tilde{\nabla }}_{\theta }\ln P_{\theta }(x_i)}$$

(10)

and we approximate $\frac{1}{N}{W}_{\theta }^f(x_i)=\frac{1}{N}w(q(x_i))$ by ${\widehat{w}}_i=\frac{1}{N}w(\frac{\mathrm{rk}(x_i)+1/2}{N})$, where rk(x_i) = |{j, f(x_j) < f(x_i)}|: it can be proven (see [1]) that ${\lim }_{N\to \infty }N{\widehat{w}}_i=W_f^{{\theta }^t}(x_i)$ (here again, we are assuming that there are no ties).

We now have an algorithm that can be used in practice if the Fisher information matrix is known.

Definition 1. The IGO update associated with parametrization θ, sample size N, step size δt and selection scheme w is given by the following update rule:

$${\theta }^{t+\delta t}={\theta }^t+\delta t{I}^{-1}({\theta }^t){\displaystyle \sum _{i=1}^N{\widehat{w}}_i\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }}.$$

(11)

We call IGO speed the vector$I^{-1}({\theta }^t){\displaystyle {\sum }_{i=1}^N{\widehat{w}}_i\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }}.$

Notice that one could start directly with the ${\widehat{w}}_i$ rather than w, as we will do later.

Replacing f by its expected value under a probability distribution P_θ and optimizing over θ using the natural gradient have already been discussed. For example, in the case of a function f defined on {0, 1}ⁿ, IGO with the Bernoulli distributions yields the algorithm, PBIL[9]. Another similar approach (stochastic relaxation) is given in [10]. For a continuous function, as we will see later, the IGO framework recovers several known ranked-based natural gradient algorithms, such as pure rank-μ CMA-ES [11], xNES or SNES (Separable Natural Evolution Strategies) [12]. See [13] or [14] for other, not necessarily gradient-based, optimization algorithms on manifolds.

2.3. Geodesic IGO

Although the IGO flow associated with a family of probability distributions is intrinsic (it only depends on the family itself, not the parametrization we choose for it), the IGO update is not. However, the difference between two steps of IGO that differ only by the parametrization is only O(δt²), whereas the different between two vanilla gradient descents with different parametrizations is O(δt).

Intuitively, the reason for this difference is that two IGO algorithms start at the same point and follow “straight lines” with the same initial speed, but the definition of “straight lines” changes with the parametrization.

For instance, in the case of Gaussian distributions, let us consider two different IGO updates with Gaussian distributions in dimension one, the first one with parametrization (μ, σ) and the second one with parametrization (μ, c := σ²). We suppose that the IGO speed for the first algorithm is $(\dot{\mu },\dot{\sigma })$. The corresponding IGO speed in the second parametrization is given by the identity $\dot{c}=2\sigma \dot{\sigma }$. Therefore, the first algorithm gives the standard deviation ${\sigma }_{\mathrm{new},1}={\sigma }_{\mathrm{old}}+\delta t\dot{\sigma }$ and the variance $c_{\mathrm{new},1}={({\sigma }_{\mathrm{new},1})}^2=c_{\mathrm{old}}+2\delta t{\sigma }_{\mathrm{old}}\dot{\sigma }+\delta t^2{\dot{\sigma }}^2=c_{\mathrm{new},2}+\delta t^2{\dot{\sigma }}^2$.

The geodesics of a Riemannian manifold are the generalization of the notion of a straight line: they are curves that locally minimize length. In particular, given two points a and b on the Riemannian manifold M, the shortest path from a to b is always a geodesic (the converse is not true, though). The notion will be explained precisely in Section 3, but let us define the geodesic IGO algorithm, which follows the geodesics of the manifold instead of following the straight lines for an arbitrary parametrization.

Definition 2 (GIGO). The geodesic IGO update (GIGO) associated with sample size N, step size δt and selection scheme w is given by the following update rule:

$${\theta }^{t+\delta t}={\exp }_{{\theta }^t}(Y\delta t)$$

(12)

where:

$$Y=I^{-1}({\theta }^t){\displaystyle \sum _{i=1}^N{\widehat{w}}_i\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }},$$

(13)

is the IGO speed and${\exp }_{{\theta }^t}$ is the exponential of the Riemannian manifold Θ. Namely, ${\exp }_{{\theta }^t}(Y\delta t)$ is the endpoint of the geodesic of Θ starting at θ^t, with initial speed Y, after a time δt. By definition, this update does not depend on the parametrization θ.

Notice that while the GIGO update is compatible with the IGO flow (in the sense that when δt → 0 and N → ∞, a parameter θ^t updated according to the GIGO algorithm is a solution of Equation (9), the equation defining the IGO flow), it not necessarily an IGO update. More precisely, the GIGO update is an IGO update if and only if the geodesics of Θ are straight lines for some parametrization (by Beltrami’s theorem, this is equivalent to Θ having constant curvature).

This is a particular case of a retraction [14]: a map from the tangent bundle of a manifold to the manifold itself satisfying a rigidity condition. Arguably, the Riemannian exponential is the most natural retraction, since it depends only on the Riemannian manifold itself and not on any decomposition. However, in general, the geodesics are difficult to compute.

In the next section, we will state Noether’s theorem, which will be our main tool to compute the GIGO update for Gaussian distributions.

3. Riemannian Geometry, Noether’s Theorem

3.1. Riemannian Geometry

The goal of this section is to state Noether’s theorem. See [15] for the proofs and [16] or [17] for a more detailed presentation. Noether’s theorem states that if a system has symmetries, then there are invariants attached to these symmetries. Firstly, we need some definitions.

Definition 3 (Motion in a Lagrangian system). Let M be a differentiable manifold, TM the set of tangent vectors on M (a tangent vector is identified by the point at which it is tangent and a vector in the tangent space) and$\begin{array}{r}\mathcal{L}:TM\to ℝ\\ (q,v)\mapsto \mathcal{L}(q,v)\end{array}$ differentiable function called the Lagrangian function (in general, it could depend on t). A “motion in the Lagrangian system (M, $ℒ$) from x to y” is map γ : [t₀, t₁] → M, such that:

• γ(t₀) = x
• γ(t₁) = y
• γ is a local extremum of the functional:

  $$\mathrm{\Phi }(\gamma )={\displaystyle {\int }_{t_0}^{t_1}ℒ(\gamma (t),\dot{\gamma }(t))\mathrm{d}t,}$$

  (14)

  among all curves c: [t₀, t₁] → M, such that c(t₀) = x, and c(t₁) = y.

For example, when (M, g) is a Riemannian manifold, the length of a curve γ between γ(t₀) and γ(t₁) is:

$${\displaystyle {\int }_{t_0}^{t_1}\sqrt{g(\dot{\gamma }(t),\dot{\gamma }(t))}\mathrm{d}t.}$$

(15)

The curves that follow the shortest path between two points x, y ∈ M are therefore the minima γ of the functional (15), such that γ(t₀) = x and γ(t₁) = y, and the corresponding Lagrangian function is $(q,v)\mapsto \sqrt{g(v,v)}$. However, any curve following the shortest trajectory will have minimum length. For example, if γ₁: [a, b] → M is a curve of the shortest path, so is γ₂: t ↦ γ₁(t²): these two curves define the same trajectory in M, but they do not travel along this trajectory at the same speed. This leads us to the following definition:

Definition 4 (Geodesics). Let I be an interval of ℝ and (M, g) be a Riemannian manifold. A curve γ: I → M is called a geodesic if for all$t_0,t_1\in I,\gamma {|}_{[t_0,t_1]}$ is a motion in the Lagrangian system (M, $ℒ$) from γ(t₀) to γ(t₁), where:

$$ℒ(\gamma )={\displaystyle {\int }_{t_0}^{t_1}g(\dot{\gamma }(t),\dot{\gamma }(t))}\mathrm{d}t$$

(16)

It can be shown (see [16]) that geodesics are curves that locally minimize length, with constant velocity, in the sense that $\frac{\mathrm{d}g(\dot{\gamma }(t),\dot{\gamma }(t))}{\mathrm{d}t}=0$. In particular, given a starting point and a starting speed, the geodesic is unique. This motivates the definition of the exponential of a Riemannian manifold.

Definition 5. Let (M, g) be a Riemannian manifold. We call the exponential of M the application:

$$\begin{array}{c}\exp :TM\to M\\ (x,v)\mapsto {\exp }_x(v),\end{array}$$

such that for any x ∈ M, if γ is the geodesic of M satisfying γ(0) = x and γ′(0) = v, then exp_x(v) = γ(1).

In order to find an extremal of a functional, the most commonly-used result is called the “Euler–Lagrange equations” (see [15], for example); a motion γ in the Lagrangian system (M, $ℒ$) must satisfy:

$$\frac{\partial ℒ}{\partial x}(\gamma (t))-\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial ℒ}{\partial \dot{x}}(\dot{\gamma }(t))\right)=0.$$

(17)

By applying this equation with the Lagrangian given by (16), it is possible to show that the geodesics of a Riemannian manifold follow the “geodesic equations”:

$${\ddot{x}}^k+{\mathrm{\Gamma }}_{ij}^k{\dot{x}}^i{\dot{x}}^j=0,$$

(18)

where the

$${\mathrm{\Gamma }}_{ij}^k=\frac{1}{2}glk\left(\frac{\partial g_{jl}}{\partial q_i}+\frac{\partial g_{li}}{\partial q_j}-\frac{\partial g_{ij}}{\partial q_l}\right)$$

(19)

are called “Christoffel symbols” of the metric g. However, these coefficients are tedious (and sometimes difficult) to compute, and (18) is a second order differential equation. Noether’s theorem will give us a first order equation to compute the geodesics.

3.2. Noether’s Theorem

Definition 6. Let h: M → M, a diffeomorphism. We say that the Lagrangian system (M, $ℒ$) admits the symmetry h if for any (q, v) ∈ TM,

$$ℒ(h(q),\mathrm{d}h(v))=ℒ(q,v),$$

(20)

where dh is the differential of h.

If M is clear in the context, we will sometimes say that $ℒ$ is invariant under h.

An example will be given in the proof of Theorem 3.

We can now state Noether’s theorem (see, for example, [15]).

Theorem 1 (Noether’s Theorem). If the Lagrangian system (M, $ℒ$) admits the one-parameter group of symmetries h^s: M → M, s ∈ ℝ, then the following quantity remains constant during motions in the system (M, $ℒ$). Namely,

$$I(\gamma (t),\dot{\gamma }(t))=\frac{\partial ℒ}{\partial v}\left(\frac{\mathrm{d}{h}^s(\gamma (t))}{\mathrm{d}s}{|}_{s=0}\right)$$

(21)

does not depend on t if γ is a motion in (M, $ℒ$).

Now, we are going to apply this theorem to our problem: computing the geodesics of Riemannian manifolds of Gaussian distributions.

4. GIGO in ${\tilde{\mathbb{G}}}_d$

If we force the covariance matrix to be either diagonal or proportional to the identity matrix, the geodesics have a simple expression that we give below. In the former case, the manifold we are considering is ${({\mathbb{G}}_1)}^d$, and in the latter case, it is ${\tilde{\mathbb{G}}}_d$.

The geodesics of ${({\mathbb{G}}_1)}^d$ are given by:

Proposition 2. Let M be a Riemannian manifold; let d ∈ ℕ; let Φ be the Riemannian exponential of M^d; and let φ be the Riemannian exponential of M. We have:

$${\mathrm{\Phi }}_{(x_1,\dots ,x_n)}((v_1,\dots ,v_n))=({\varphi }_{x_1}(v_1),\dots ,{\varphi }_{xn}(v_n))$$

(22)

In particular, knowing the geodesics of${\mathbb{G}}_1$ is enough to compute the geodesics of${({\mathbb{G}}_1)}^d$.

This is true, because a block of the product metric does not depend on variables of the other blocks. Consequently, a GIGO update with a diagonal covariance matrix with the sample (x_i) is equivalent to d separate one-dimensional GIGO updates using the same samples. Moreover, ${\mathbb{G}}_1\cong {\tilde{\mathbb{G}}}_1$, the geodesics of which are given below.

We will show that ${\tilde{\mathbb{G}}}_d$ and the “hyperbolic space”, of which the geodesics are known, are isometric.

4.1. Preliminaries: Poincaré Half-Plane, Hyperbolic Space

In dimension two, the hyperbolic space is called the “hyperbolic plane” or the Poincaré half-plane. We recall its definition:

Definition 7 (Poincaré half-plane). We call the “Poincaré half-plane” the Riemannian manifold:

$$\mathscr{H}=\{(x,y)\in {ℝ}^2,y>0\},$$

with the metric$d{s}^2=\frac{d{x}^2+d{y}^2}{y^2}$.

We also recall the expression of its geodesics (see, for example, [18]):

Proposition 3 (Geodesics of the Poincaré half-plane). The geodesics of the Poincaré half-plane are exactly the:

$$t\mapsto (\mathrm{Re}(z(t)),\mathrm{Im}(z(t))),$$

where:

$$z(t)=\frac{ai{e}^{vt}+b}{ci{e}^{vt}+d},$$

(23)

with ad − bc = 1 and v > 0.

The geodesics are half-circles perpendicular to the line y = 0 and vertical lines, as shown in Figure 1 below.

The generalization to the higher dimension is the following:

Definition 8 (Hyperbolic space). We call the “hyperbolic space of dimension n” the Riemannian manifold:

$${\mathscr{H}}_n=\{(x_1,\dots ,x_{n-1,}y)\in {ℝ}^n,y>0\},$$

with the metric$d{s}^2=\frac{d{x}_1^2+\dots +d{x}_{n-1}^2+d{y}^2}{y^2}$ (or equivalently, the metric given by the matrix$\mathrm{Diag}\left(\frac{1}{y^2}\right)$).

The Lagrangian for the geodesics is invariant under all translations along the x_i, so by Noether’s theorem, its geodesics stay in a plane containing the direction y and the initial speed. The induced metric on this plane is the metric of the Poincaré half-plane. The geodesics are therefore given by the following proposition:

Proposition 4 (Geodesics of the hyperbolic space). If γ : t ⟼ (x₁(t), …, x_n−1(t), y(t)) = (x(t), y(t)) is a geodesic of $\mathscr{H}$_n, then there exists a, b, c, d ∈ ℝ, such that ad − bc = 1, and v > 0, such that$x(t)=x(0)+\frac{{\dot{x}}_0}{\Vert {\dot{x}}_0\Vert }\tilde{x}(t),y(t)=\mathrm{Im}({\gamma }_{ℂ}(t))$, with$\tilde{x}(t)=\mathrm{Re}({\gamma }_{ℂ}(t))$ and:

$${\gamma }_{ℂ}(t):=\frac{ai{e}^{vt}+b}{ci{e}^{vt}+d}.$$

(24)

4.2. Computing the GIGO Update in ${\tilde{\mathbb{G}}}_d$

If we want to implement the GIGO algorithm in ${\tilde{\mathbb{G}}}_d$, we need to compute the natural gradient in ${\tilde{\mathbb{G}}}_d$ and to be able to compute the Riemannian exponential of ${\tilde{\mathbb{G}}}_d$.

Using Proposition 1, we can compute the metric of ${\tilde{\mathbb{G}}}_d$ in the parametrization (μ, σ) ⟼ N(μ, σ²I). We find:

$$\left(\begin{array}{cccc}\frac{1}{{\sigma }^2}& 0& \dots & 0\\ 0& \ddots & \ddots & ⋮\\ ⋮& \ddots & \frac{1}{{\sigma }^2}& 0\\ 0& \dots & 0& \frac{2d}{{\sigma }^2}\end{array}\right).$$

(25)

Since this matrix is diagonal, it is easy to invert, and we immediately have the natural gradient and, consequently, the IGO speed.

Proposition 5. In${\tilde{\mathbb{G}}}_d$, the IGO speed Y is given by:

$$Y_{\mu }={\displaystyle \sum _i{\widehat{w}}_i(x_i-\mu ),}$$

(26)

$$Y_{\sigma }={\displaystyle \sum _i{\widehat{w}}_i\left(\frac{{(x_i-\mu )}^T(x_i-\mu )}{2d\sigma }-\frac{\sigma }{2}\right).}$$

(27)

Proof. We recall the IGO speed is defined by $Y=I^{-1}({\theta }^t){\displaystyle {\sum }_{i=1}^N{\widehat{w}}_i\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }}$. Since $P_{\mu ,\sigma }(x)={(2\pi {\sigma }^2)}^{-d/2}\exp (-\frac{{(x-\mu )}^T(x-\mu )}{2{\sigma }^2})$, we have:

$$\begin{array}{c}\frac{\partial \ln P_{\mu ,\sigma }(x)}{\partial \mu }=x-\mu ,\\ \frac{\partial \ln P_{\mu ,\sigma }(x)}{\partial \sigma }=-\frac{d}{\sigma }+\frac{{(x-\mu )}^T(x-\mu )}{{\sigma }^3}.\end{array}$$

The result follows. □

The metric defined by Equation (25) is not exactly the metric of the hyperbolic space, but with the substitution $\mu \leftarrow \frac{\mu }{\sqrt{2d}}$, the metric becomes $\frac{2d}{{\sigma }^2}I$, which is proportional to the metric of the hyperbolic space and, therefore, defines the same geodesics.

Theorem 2 (Geodesics of ${\tilde{\mathbb{G}}}_d$). If$\gamma :t\mapsto \mathcal{N}(\mu (t),\sigma {(t)}^{2}I)$ is a geodesic of${\tilde{\mathbb{G}}}_d$, then there exists a, b, c, d ∈ ℝ, such that ad − bc = 1, and v > 0, such that:$\mu (t)=\mu (0)+\sqrt{2d}\frac{{\dot{\mu }}_0}{\Vert {\dot{\mu }}_0\Vert }\tilde{r},\sigma (t)=\mathrm{Im}({\gamma }_{ℂ}(t))$, with$\tilde{r}(t)=\mathrm{Re}({\gamma }_{ℂ}(t))$and

$${\gamma }_{ℂ}(t):=\frac{ai{e}^{vt}+b}{ci{e}^{vt}+d}.$$

(28)

Now, in order to implement the corresponding GIGO algorithm, we only need to be able to find the coefficients a, b, c, d, v corresponding to an initial position (μ₀, σ₀) and an initial speed $({\dot{\mu }}_0,{\dot{\sigma }}_0)$. This is a tedious but easy computation, the result of which is given in Proposition 17.

The pseudocode of GIGO in ${\tilde{\mathbb{G}}}_d$ is also given in the Appendix: it is obtained by concatenating

Algorithm 1. For all algorithms.

and

Algorithm 7. GIGO in ${\tilde{\mathbb{G}}}_d$, one step.

(Proposition 17 and the pseudocode in the Appendix allow the metric to be slightly modified; see Section 6.2).

5. GIGO in ${\tilde{\mathbb{G}}}_d$

5.1. Obtaining a First Order Differential Equation for the Geodesics of ${\mathbb{G}}_d$

In the case where both the covariance matrix and the mean can vary freely, the equations of the geodesics have been computed in [4] and [5]. However, these articles start with the equations of the geodesics obtained with the Christoffel symbols, then partially integrate them. These equations are in fact a consequence of Noether’s theorem and can be found directly.

Theorem 3. Let$\gamma :t\mapsto \mathcal{N}({\mu }_t,{\sum }_t)$ be a geodesic of${\mathbb{G}}_d$. Then, the following quantities do not depend on t:

$$J_{\mu }={\displaystyle {\sum }_t^{-1}{\dot{\mu }}_t,}$$

(29)

$$J_{\sum }={\displaystyle {\sum }_t^{-1}({\dot{\mu }}_t{\mu }_t^T+{\dot{\sum }}_t).}$$

(30)

Proof. This is a direct application of Noether’s theorem, with suitable groups of diffeomorphisms. By Proposition 1, the Lagrangian associated with the geodesics of ${\mathbb{G}}_d$ is:

$$\mathcal{L}(\mu ,\sum ,\dot{\mu },\dot{\sum })={\dot{\mu }}^T{\sum }^{-1}\dot{\mu }+\frac{1}{2}\mathrm{tr}(\dot{\sum }{\sum }^{-1}\dot{\sum }{\sum }^{-1}).$$

(31)

Its derivative is:

$$\frac{\partial ℒ}{\partial \dot{\theta }}=[(h,H)\mapsto 2{\dot{\mu }}^T{\sum }^{-1}h+\mathrm{tr}(H{\sum }^{-1}\dot{\sum }{\sum }^{-1})]$$

(32)

Let us show that this Lagrangian is invariant under affine changes of basis (thus illustrating Definition 6).

The general form of an affine change of basis is ϕ_μ0,A : (μ, Σ) ⟼ (Aμ + μ₀, AΣA^T ), with μ₀ ∈ ℝ^d and A ∈ GL_d(ℝ).

We have:

$$\mathcal{L}({\varphi }_{\mu 0,A}(\mu ,\sum ),d{\varphi }_{\mu 0,A}(\dot{\mu },\dot{\sum }))={\stackrel{.}{\overline{A\mu }}}^T{(A\sum A^T)}^{-1}\stackrel{.}{\overline{A\mu }}+\frac{1}{2}\mathrm{tr}(\overline{A{\displaystyle \sum }{A}^T}{(A\sum A^T)}^{-1}\overline{A{\displaystyle \sum }{A}^T}{(A\sum A^T)}^{-1}.)$$

(33)

and since $\stackrel{.}{\overline{A\mu }}=A\dot{\mu }$ and, $\overline{A{\displaystyle \sum }{A}^T}(A\sum A^T)$ we find easily that:

$$ℒ({\varphi }_{\mu 0,A}(\mu ,\mathrm{\Sigma }),d{\varphi }_{\mu 0,A}(\dot{\mu },\dot{\mathrm{\Sigma }}))=ℒ(\mu ,\mathrm{\Sigma },\dot{\mu },\dot{\mathrm{\Sigma }}),$$

(34)

or in other words: $ℒ$ is invariant under ϕ_μ0,A for any μ₀ ∈ ℝ^d, A ∈ GL_d(ℝ).

In order to use Noether’s theorem, we also need one-parameter groups of transformations. We choose the following:

• Translations of the mean vector. For any i ∈ [1, d], $\mathrm{let}{h}_i^s:(\mu ,\mathrm{\Sigma })\mapsto (\mu +s{e}_i,\mathrm{\Sigma })$, where e_i is the i-th basis vector. We have $\frac{\mathrm{d}{h}_i^s}{\mathrm{d}s}{|}_{s=0}=(e_i,0)$, so by Noether’s theorem,

  $$\frac{\partial ℒ}{\partial \dot{\theta }}(e_{i,}0)=2{\dot{\mu }}^T{\mathrm{\Sigma }}^{-1}{e}_i=2e_i^T{\mathrm{\Sigma }}^{-1}\dot{\mu }$$

  remains constant for all i. The fact that J_μ is an invariant immediately follows.
• Linear base changes. For any i, j ∈ [1, d], $\mathrm{let}{h}_{i,j}^s:(\mu ,\mathrm{\Sigma })\mapsto (\exp (s{E}_{ij})\mu ,\exp (s{E}_{ij})\mathrm{\Sigma }\exp (s{E}_{ij})),$, where E_ij is the matrix with a one at position (i, j) and zeros elsewhere. We have:

  $$\frac{\mathrm{d}{h}_{E_{ij}}^s}{\mathrm{d}s}{|}_{s=0}=(E_{ij\mu },E_{ij}\mathrm{\Sigma }+E_{ji}\mathrm{\Sigma }).$$

Therefore, by Noether’s theorem, we then obtain the following invariants:

$$J_{ij}:=\frac{\partial ℒ}{\partial \dot{\theta }}(E_{ij\mu },E_{ij}\mathrm{\Sigma }+E{\mathrm{\Sigma }}_{ij})$$

(35)

$$=2{\dot{\mu }}^T{\mathrm{\Sigma }}^{-1}{E}_{ij\mu }+\mathrm{tr}((E_{ij\mu }\mathrm{\Sigma }+\mathrm{\Sigma }{E}_{ji}){\mathrm{\Sigma }}^{-1}\dot{\mathrm{\Sigma }}{\mathrm{\Sigma }}^{-1})$$

(36)

$$=2{({\mathrm{\Sigma }}^{-1}\dot{\mu })}^T{E}_{ij}\mu +\mathrm{tr}(E_{ij}\dot{\mathrm{\Sigma }}{\mathrm{\Sigma }}^{-1})+\mathrm{tr}(E_{ji}{\mathrm{\Sigma }}^{-1}\dot{\mathrm{\Sigma }})$$

(37)

$$=2{(J_{\mu }{\mu }^T)}_{ij}+2{({\mathrm{\Sigma }}^{-1}\dot{\mathrm{\Sigma }})}_{ij},$$

(38)

and the coefficients of J_Σ in (30) are the (J_ij/2).

This leads us to first order equations satisfied by the geodesics mentioned in [4–6].

Theorem 4 (GIGO-Σ). $t7\mapsto \mathcal{N}({\mu }_t,{\mathrm{\Sigma }}_t)$ is a geodesic of${\mathbb{G}}_d$ if and only if μ : t 7↦μ_t and Σ : t ↦Σ_t satisfy the equations:

$${\dot{\mu }}_t={\mathrm{\Sigma }}_t{J}_{\mu }$$

(39)

$${\dot{\mathrm{\Sigma }}}_t={\mathrm{\Sigma }}_t(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }_t^T)={\mathrm{\Sigma }}_t{J}_{\mathrm{\Sigma }}-{\dot{\mu }}_t{\mu }_t^T,$$

(40)

where:

$$J_{\mu }={\displaystyle \sum {}_0^{-1}{\dot{\mu }}_0,}$$

and:

$$J_{\mathrm{\Sigma }}={\mathrm{\Sigma }}_0^{-1}\left({\dot{\mu }}_0{\mu }_0^T+{\dot{\mathrm{\Sigma }}}_0\right).$$

Proof. This is an immediate consequence of Proposition 3.

These equations can be solved analytically (see [5]); however, usually, that is not the case, and they have to be solved numerically, for example with the Euler method (the corresponding algorithm, which we call GIGO-Σ, is described in the Appendix). The goal of the remainder of the subsection is to show that having to use the Euler method is fine.

To avoid confusion, we will call the step size of the GIGO algorithm (δt in Proposition 2) “GIGO step size” and the step size of the Euler method (inside a step of the GIGO algorithm) “Euler step size”.

Having to solve our equations numerically brings two problems:

The first one is a theoretical problem: the main reason to study GIGO is its invariance under reparametrization of θ, and we lose this invariance property when we use the Euler method. However, GIGO can get arbitrarily close to invariance by decreasing the Euler step size. In other words, the difference between two different IGO algorithms is O(δt²), and the difference between two different implementations of the GIGO algorithm is O(h²), where h is the Euler step size; it is easier to reduce the latter. Still, without a closed form for the geodesics of ${\mathbb{G}}_d$, the GIGO update is rather expensive to compute, but it can be argued that most of the computation time will still be the computation of the objective function f.

The second problem is purely numerical: we cannot guarantee that the covariance matrix remains positive-definite along the Euler method. Here, apart from finding a closed form for the geodesics, we have two solutions.

We can enforce this a posteriori: if the covariance matrix we find is not positive-definite after a GIGO step, we repeat the failed GIGO step with a reduced Euler step size (in our implementation, we divided it by four; see

Algorithm 2. GIGO Update, one step, updating the covariance matrix.

in the Appendix.).

The other solution is to obtain differential equations on a square root of the covariance matrix (any matrix A, such that Σ = AA^T ).

Theorem 5 (GIGO-A). If μ : t ↦μ_t and A : t ↦A_t satisfy the equations:

$${\dot{\mu }}_t=A_t{A}_t^T{J}_{\mu },$$

(41)

$${\dot{A}}_t=\frac{1}{2}{\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }_t^T\right)}^T{A}_t,$$

(42)

where:

$$J_{\mu }={(A_0^{-1})}^T{A}_0^{-1}{\mu }_0$$

and:

$$J_{\mathrm{\Sigma }}={(A_0^{-1})}^T{A}_0^{-1}({\dot{\mu }}_0{\mu }_t^T+{\dot{A}}_0{A}_0^T+A_0{\dot{A}}_0^T),$$

then$t\mapsto \mathcal{N}({\mu }_t,A_t{A}_t^T)$ is a geodesic of${\mathbb{G}}_d$.

Proof. This is a simple rewriting of Theorem 4: if we write Σ := AA^T, we find that J_μ and J_Σ are the same as in Theorem 4, and we have:

$$\dot{\mu }=\mathrm{\Sigma }{J}_{\mu },$$

and:

$$\begin{array}{c}\dot{\mathrm{\Sigma }}=(\dot{A}{A}^T+A{\dot{A}}^T)=\frac{1}{2}{\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)}^T+A{A}^T+\frac{1}{2}A{A}^T\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)\\ =\frac{1}{2}{\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)}^T\mathrm{\Sigma }+\frac{1}{2}\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)=\frac{1}{2}\mathrm{\Sigma }\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)+\frac{1}{2}{\left[\mathrm{\Sigma }\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)\right]}^T.\end{array}$$

By Theorem 4, Σ(J_Σ − J_μμ^T) is symmetric (since $\dot{\mathrm{\Sigma }}$ has to be symmetric). Therefore, we have $\dot{\mathrm{\Sigma }}\left(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T\right)$, and the result follows. □

Notice that Theorem 5 gives an equivalence, whereas Theorem 4 does not. The reason is that the square root of a symmetric positive-definite matrix is not unique. Still, it is canonical; see the discussion in Section 6.1.2.

As for Theorem 4, we can solve Equations (41) and (42) numerically, and we obtain another algorithm (

Algorithm 3. GIGO Update, one step, updating a square root of the covariance matrix.

in the Appendix; we will call it GIGO-A), with a behavior similar to the previous one (with Equations (39) and (40)). For both of them, numerical problems can arise when the covariance matrix is almost singular.

We have not managed to find any example where one of these two algorithms converged to the minimum of the objective function, whereas the other did not, and their behavior is almost the same.

More interestingly, the performances of these two algorithms are also the same as the performances of the exact GIGO algorithm, using the equations of Section 5.2.

Notice that even though GIGO-A directly maintains a square root of the covariance matrix, which makes sampling new points easier (to sample a point from $\mathcal{N}\left(\mu ,\mathrm{\Sigma }\right)$, a square root of Σ is needed), both GIGO-Σ and GIGO-A still have to invert the covariance matrix (or its square root) at each step, which is as costly as the decomposition, so one of these algorithms is roughly as expensive to compute as the other.

5.2. Explicit Form of the Geodesics of ${\mathbb{G}}_d$ (from [5])

We now give the exact geodesics of ${\mathbb{G}}_d$: the following results are a rewriting of Theorem 3.1 and its first corollary in [5].

Theorem 6. Let$\left({\dot{\mu }}_0,{\dot{\mathrm{\Sigma }}}_0\right)\in T_{\mathcal{N}}\left(0,I\right){\mathbb{G}}_d$. The geodesic of${\mathbb{G}}_d$ starting from$\mathcal{N}(0,1)$ with initial speed$\left({\dot{\mu }}_0,{\dot{\mathrm{\Sigma }}}_0\right)$ is given by:

$${\exp }_{\mathcal{N}(0,I)}(s{\dot{\mu }}_0,s{\dot{\mathrm{\Sigma }}}_0)=\mathcal{N}\left(2R(s)\mathrm{sh}\left(\frac{s\mathbb{G}}{2}\right){\mathbb{G}}^{-}{\dot{\mu }}_0,R(s)R{(s)}^T\right),$$

(43)

where exp is the Riemannian exponential of${\mathbb{G}}_d$, G is any matrix satisfying:

$${\mathbb{G}}^2={\dot{\mathrm{\Sigma }}}_0^2+2{\dot{\mu }}_0{\dot{\mu }}_0^T,$$

(44)

$$R(s)={\left({\left(\mathrm{ch}(\frac{s\mathbb{G}}{2})-{\dot{\mathrm{\Sigma }}}_0{\mathbb{G}}^{-}\mathrm{sh}(\frac{s\mathbb{G}}{2})\right)}^{-1}\right)}^T$$

(45)

and G⁻ is a pseudo-inverse of G

In [5], the existence of G (as a square root of ${\dot{\mathrm{\Sigma }}}_0^2+2{\dot{\mu }}_0{\dot{\mu }}_0^T$) is proven. Notice that, anyway, in the expansions of (43) and (45), only even powers of G appear.

Additionally, since, for all A ∈ GL_d(ℝ), for all μ₀ ∈ ℝ^d, the application:

$$\begin{array}{l}\varphi :{\mathbb{G}}_d\to {\mathbb{G}}_d\\ \mathcal{N}(\mu ,\mathrm{\Sigma })\mapsto \mathcal{N}(A\mu +{\mu }_0,A\mathrm{\Sigma }{A}^T)\end{array}$$

(46)

preserves the geodesics, we find the general expression for the geodesics of ${\mathbb{G}}_d$.

Corollary 1. Let μ₀ ∈ ℝ^d, A ∈ GL_d(ℝ) and$({\dot{\mu }}_0,{\dot{\mathrm{\Sigma }}}_0)\in T_{\mathcal{N}({\mu }_{0.},A_0{A}_0^T)}{\mathbb{G}}_d$. The geodesic of${\mathbb{G}}_d$ starting from$\mathcal{N}\left(\mu ,\mathrm{\Sigma }\right)$ with initial speed$({\dot{\mu }}_0,{\dot{\mathrm{\Sigma }}}_0)$ is given by:

$${\exp }_{\mathcal{N}({\mu }_0,A_0{A}_0^T)}(s{\dot{\mu }}_0,s{\dot{\mathrm{\Sigma }}}_0)=\mathcal{N}({\mu }_1,A_1{A}_1^T),$$

(47)

$${\mu }_1=2A_{0}R(s)\mathrm{sh}(\frac{s\mathbb{G}}{2}){\mathbb{G}}^{-}{A}_0^{-1}{\dot{\mu }}_0+{\mu }_0,$$

(48)

$$A_1=A_{0}R(s),$$

(49)

where exp is the Riemannian exponential of${\mathbb{G}}_d$, G is any matrix satisfying:

$${\mathbb{G}}^2=A_0^{-1}\left({\dot{\mathrm{\Sigma }}}_0{\mathrm{\Sigma }}_0^{-1}{\dot{\mathrm{\Sigma }}}_0+2{\dot{\mu }}_0{\dot{\mu }}_0^T\right){(A_0^{-1})}^T,$$

(50)

$$R(s)={\left({\left(\mathrm{ch}(\frac{s\mathbb{G}}{2})-A_0^{-1}{\dot{\mathrm{\Sigma }}}_0{(A_0^{-1})}^T{\mathbb{G}}^{-}\mathrm{sh}(\frac{s\mathbb{G}}{2})\right)}^{-1}\right)}^T,$$

(51)

and G⁻ is a pseudo-inverse of G.

It should be noted that the final values for mean and covariance do not depend on the choice of G as a square root of:

$$A_0^{-1}\left({\dot{\mathrm{\Sigma }}}_0{\mathrm{\Sigma }}_0^{-1}{\dot{\mathrm{\Sigma }}}_0+2{\dot{\mu }}_0{\dot{\mu }}_0^T\right){(A_0^{-1})}^T.$$

The reason for this is that ch(G) is a Taylor series in G², and so are sh(G)G⁻ and G⁻sh(G).

For our practical implementation, we actually used these Taylor series instead of the expression of the corollary

6. Comparing GIGO, xNES and Pure Rank-μ CMA-ES

6.1. Definitions

In this section, we recall the xNES and pure rank-μ CMA-ES, and we describe them in the IGO framework, thus allowing a reasonable comparison with the GIGO algorithms.

6.1.1. xNES

We recall a restriction of the xNES algorithm, introduced in [19] (this restriction is sufficient to describe the numerical experiments in [19]).

Definition 9 (xNES algorithm). The xNES algorithm with sample size N, weights w_i and learning rates η_μ and η_Σ updates the parameters μ ∈ ℝ^d, A ∈ M_d(ℝ) with the following rule: At each step, N points x₁, …, x_N are sampled from the distribution.$\mathcal{N}(\mu ,A{A}^T)$. Without loss of generality, we assume f(x₁) < … < f(x_N). The parameter is updated according to:

$$\begin{array}{c}\mu \leftarrow \mu +{\eta }_{\mu }A{\mathbb{G}}_{\mu },\\ A\leftarrow A\exp ({\eta }_{\mathrm{\Sigma }}{\mathbb{G}}_M/2),\end{array}$$

where, setting z_i = A⁻¹(x_i − μ):

$$\begin{array}{c}{\mathbb{G}}_{\mu }={\displaystyle \sum _{i=1}^N{w}_i{z}_i,}\\ {\mathbb{G}}_M={\displaystyle \sum _{i=1}^N{w}_i(z_i{z}_i^T-I).}\end{array}$$

The more general version decomposes the matrix A as σB, where det B = 1, and uses two different learning rates for σ and for B. We gave the version where these two learning rates are equal (in particular, for the default parameters in [19], these two learning rates are equal). This restriction of the xNES algorithm can be described in the IGO framework, provided all of the learning rates are equal (most of the elements of the proof can be found in [19] (the proposition below essentially states that xNES is a natural gradient update) or in [1]):

Proposition 6 (xNES as IGO). The xNES algorithm with sample size N, weights w_i and learning rates η_μ = η_Σ = δt coincides with the IGO algorithm with sample size N, weights w_i, step size δt and in which, given the current position (μ_t, A_t), the set of Gaussians is parametrized by:

$${\varphi }_{{\mu }_t,A_t}:(\delta ,M)\mapsto \mathcal{N}\left({\mu }_t+A_t\delta ,\left(A_t\exp (\frac{1}{2}M)\right){\left(A_t\exp (\frac{1}{2}M)\right)}^T\right),$$

with δ ∈ ℝ^m and M ∈ Sym(ℝ^m).

The parameters maintained by the algorithm are (μ, A), and the x_i are sampled from$\mathcal{N}(\mu ,A{A}^T)$.

Proof. Let us compute the IGO update in the parametrization ${\varphi }_{{\mu }_t,A_t}$: we have δ^t = 0, M^t = 0, and by using Proposition 1, we can see that for this parametrization, the Fisher information matrix at (0, 0) is the identity matrix. The IGO update is therefore,

$${(\delta ,M)}^{t+\delta t}={(\delta ,M)}^t+\delta t{Y}_{\delta }(\delta ,M)+\delta t{Y}_M(\delta ,M)=\delta t{Y}_{\delta }(\delta ,M)+\delta t{Y}_M(\delta ,M),$$

where:

$$Y_{\delta }(\delta ,M)={\displaystyle \sum _{i=1}^N{w}_i{\nabla }_{\delta }\ln (p(x_i|(\delta ,M))}$$

and:

$$Y_M(\delta ,M)={\displaystyle \sum _{i=1}^N{w}_i{\nabla }_M\ln (p(x_i|(\delta ,M))}.$$

Since tr(M) = log(det(exp(M))), we have:

$$\begin{array}{l}\ln P_{\delta ,M}(x)=-\frac{d}{2}\ln (2\pi )-\ln (\det A)-\frac{1}{2}\mathrm{tr}M-\frac{1}{2}\Vert \exp (-\frac{1}{2}M)A^{-1}(x-\mu -A\delta )|{|}^2,\\ \end{array}$$

and a straightforward computation yields:

$$Y_{\delta }(\delta ,M)={\displaystyle \sum _{i=1}^N{w}_i{z}_i}={\mathbb{G}}_{\mu },$$

and:

$$Y_M(\delta ,M)=\frac{1}{2}{\displaystyle \sum _{i=1}^N{w}_i(z_i{z}_i^T-I)={\mathbb{G}}_M.}$$

Therefore, the IGO update is:

$$\begin{array}{l}\delta (t+\delta t)=\delta (t)+\delta t{\mathbb{G}}_{\mu },\\ M(t+\delta t)=M(t)+\delta t{\mathbb{G}}_M,\end{array}$$

or, in terms of mean and covariance matrix:

$$\begin{array}{l}\mu (t+\delta t)=\mu (t)+\delta tA(t){\mathbb{G}}_{\mu }\\ A(t+\delta t)=A(t)\exp (\delta t{\mathbb{G}}_M/2),\end{array}$$

or:

$$\mathrm{\Sigma }(t+\delta t)=A(t)\exp (\delta t{\mathbb{G}}_M)A{(t)}^T.$$

This is the xNES update. □

6.1.2. Using a Square Root of the Covariance Matrix

Firstly, we recall that the IGO framework (on ${\mathbb{G}}_d$, for example) emphasizes the Riemannian manifold structure on ${\mathbb{G}}_d$. All of the algorithms studied here (including GIGO, which is not strictly speaking an IGO algorithm) define a trajectory in ${\mathbb{G}}_d$ (a new point for each step), and to go from a point θ to the next one (θ′), we follow some curve $\gamma :[0,\delta t]\to {\mathbb{G}}_d$, with γ(0) = θ, γ(δt) = θ′ and $\dot{\gamma }(0)$ given by the natural gradient $(\dot{\gamma }(0)={{\displaystyle {\sum }_{i=1}^N{\widehat{w}}_i{\tilde{\nabla }}_{\theta }{P}_{\theta }(x_i)\in T}}_{\theta }{\mathbb{G}}_d)$_.

To be compatible with this point of view, an algorithm giving an update rule for a square root (any matrix A such that Σ = AA^T: since we do not force A to be symmetric, the decomposition is not unique) of the covariance matrix A has to satisfy the following condition: for a given initial speed, the covariance matrix Σ^t+δt after one step must depend only on Σ^t and not on the square root A^t chosen for Σ^t.

The xNES algorithm does satisfy this condition: consider two xNES algorithms, with the same learning rates, respectively, at $(\mu ,A_1^t)$ and $(\mu ,A_2^t)$, with $A_1^t{(A_1^t)}^T=A_2^t{(A_2^t)}^T$ (i.e., they define the same Σ^t), using the same samples x_i to compute the natural gradient update, then we will have ${\displaystyle {\sum }_1^{t+\delta t}=}{\displaystyle \sum {}_2^{t+\delta t}}$. Using the definitions of Section 6.3, we have just shown that what we will call the “xNES trajectory” is well defined.

It is also important to notice that, in order to be well defined, a natural gradient algorithm updating a square root of the covariance matrix has to specify more conditions than simply following the natural gradient.

The reason for this is that the natural gradient is a vector tangent to ${\mathbb{G}}_d$: it lives in a space of dimension d(d + 3)/2 (the dimension of ${\mathbb{G}}_d$), whereas the vector (μ, A) lives in a space of dimension d(d + 1) (the dimension of ℝⁿ × GL_n(ℝ)), which is too large: there exists infinitely many applications t ↦A_t, such that a given curve $\gamma :t\mapsto \mathcal{N}({\mu }_t,{\mathrm{\Sigma }}_t)$ can be written $\gamma (t)=\mathcal{N}\left({\mu }_t{A}_t{A}_t^T\right)$. This is why Theorem 5 is simply an implication, whereas Theorem 4 is an equivalence.

More precisely, let us consider A in GL_d(ℝ) and v_A, $v_A^{\prime }$ two infinitesimal updates of A. Since Σ = AA^T, the infinitesimal update of Σ corresponding to $v_A^{\prime }$ (resp. $v_A^{\prime }$) is $v_{\mathrm{\Sigma }}=A{v}_A^T+v_A{A}^T$ (resp. $v_{\mathrm{\Sigma }}^{\prime }=A{v^{\prime }}_A^T+v_A^{\prime }{A}^T).$

It is now easy to see that v_A and $v_A{}^{\prime }$ define the same direction for Σ (i.e., $v_{\mathrm{\Sigma }}=v_{\mathrm{\Sigma }}{}^{\prime }$) if and only if AM^T + MA^T = 0, where $M=v_A-v_A{}^{\prime }$. This is equivalent to A⁻¹M antisymmetric.

For any A ∈ M_d(ℝ), let us denote by T_A the space of the matrices M, such that A⁻¹M is antisymmetric or, in other words, T_A := {u ∈ M_d(ℝ), Au^T + uA^T = 0}. Having a subspace S_A in direct sum with T_A for all A is sufficient (but not necessary) to have a well-defined update rule. Namely, consider the (linear) application:

$$\begin{array}{ccc}{\varphi }_A:{\mathrm{M}}_d\left(ℝ\right)& \hfill \to & {\mathrm{S}}_d\left(ℝ\right)\\ \hfill v_A& \hfill \mapsto & A{v}_A^T+v_A{A}^T\end{array},$$

sending an infinitesimal update of A to the corresponding update of Σ. It is not bijective, but as we have seen before, Ker ϕ_A = T_A, and therefore, if we have, for some U_A,

$${\mathrm{M}}_d\left(ℝ\right)=U_A\oplus T_A,$$

(52)

then φ_A|_UA is an isomorphism. Let v_Σ be an infinitesimal update of Σ. We choose the following update of A corresponding to v_Σ:

$$v_A:{({\varphi }_A|U_A)}^{-1}\left(v_{\mathrm{\Sigma }}\right).$$

(53)

Any U_A, such that U_A ⊕ T_A = M_d(ℝ), is a reasonable choice to pick v_A for a given v_Σ. The choice S_A = {u ∈ M_d(ℝ), Au^T − uA^T = 0} has an interesting additional property; it is the orthogonal of T_A for the norm:

$${\Vert v_A\Vert }_{\sum }^2:=\mathrm{Tr}(v_A^T{\sum }^{-1}{v}_A)=\mathrm{Tr}({(A^{-1}{v}_A)}^T{A}^{-1}{v}_A).$$

(54)

and consequently, it can be defined without referring to the parametrization, which makes it a canonical choice. To prove this, remark that T_A = {M ∈ M_d(ℝ), A⁻¹M antisymmetric} and S_A = {M ∈ M_d(ℝ), A⁻¹M symmetric} and that if M is symmetric and N is antisymmetric, then

$$\text{Tr}(M^{T}N)={\displaystyle \sum _{i,j=1}^d{m}_{ij}{n}_{ij}={\displaystyle \sum _{i=1}^d{m}_{ii}{n}_{ii}+{\displaystyle \sum _{1\le i\le j\le d}{m}_{ij}(n_{ij}+n_{ji})}=0.}}$$

(55)

Let us now show that this is the choice made by xNES and GIGO-A (which are well-defined algorithms updating a square root of the covariance matrix).

Proposition 7. Let A ∈ M_n(ℝ). The v_A given by the xNES and GIGO-A algorithms lies in S_A = {u ∈ M_d(ℝ), Au^T − uA^T = 0} = S_A.

Proof. For xNES, let us write $\dot{\gamma }(0)=({\upsilon }_{\mu },{\upsilon }_{\Sigma })$ and ${\upsilon }_A:=\frac{1}{2}A{G}_M$. We have $A^{-1}{\upsilon }_A=\frac{1}{2}{G}_M$, and therefore, forcing M (and G_M) to be symmetric in xNES is equivalent to A⁻¹ υ_A = (A⁻¹ υ_A)^T, which can be rewritten as $A{\upsilon }_A^T={\upsilon }_A{A}^T$. For GIGO-A, Equation (40) shows that ${\Sigma }_t(J_{\Sigma }-J_{\mu }{\mu }_t^T)$ is symmetric, and with this fact in mind, Equation (42) shows that we have $A{\upsilon }_A^T={\upsilon }_A{A}^T({\upsilon }_A\text{is}{\dot{A}}_t)$. □

6.1.3. Pure Rank-μ CMA-ES

We now recall the pure rank-μ CMA-ES algorithm. The general CMA-ES algorithm is described in [21].

Definition 10 (Pure rank-μ CMA-ES algorithm). The pure rank-μ CMA-ES algorithm with sample size N, weights w_i and learning rates η_μ and η_Σ is defined by the following update rule: At each step, N points x₁, …, x_N are sampled from the distribution$N(\mu ,\Sigma )$. Without loss of generality, we assume f)x₁) < … < f(x_N). The parameter is updated according to:

$$\mu \leftarrow \mu +{\eta }_{\mu }{\displaystyle \sum _{i=1}^N{w}_i(x_i-\mu )},$$

$$\Sigma \leftarrow \Sigma +{\eta }_{\Sigma }{\displaystyle \sum _{i=1}^N{w}_i((x_i-\mu ){(x_i-\mu )}^T-\Sigma )}.$$

The pure rank-μ CMA-ES can also be described in the IGO framework; see, for example, [20].

Proposition 8 (Pure rank-μ CMA-ES as IGO). The pure rank-μ CMA-ES algorithm with sample size N, weights w_i and learning rates η_μ = η_Σ = δt coincides with the IGO algorithm with sample size N, weights w_i, step size δt and the parametrization (μ, Σ).

6.2. Twisting the Metric

As we can see, the IGO framework does not allow one to recover the learning rates for xNES and pure rank-μ CMA-ES, which is a problem, since usually, the covariance learning rate is set much smaller than the mean learning rate (see either [19] or [21]).

A way to recover these learning rates is to incorporate them directly into the metric (see also blockwise GIGO, in Section 6.4). More precisely:

Definition 11 (Twisted Fisher metric). Let η_μ, η_Σ ∈ ℝ, and let (P_θ)_θ∈Θ be a family of normal probability distributions: P_θ = N (μ(θ), Σ(θ)), with μ and Σ C¹. We call the “(η_μ, η_Σ)-twisted Fisher metric” the metric defined by:

$$I_{i,j}({\eta }_{\mu },{\eta }_{\mathrm{\Sigma }})(\theta )=\frac{1}{{\eta }_{\mu }}\frac{\partial {\mu }^T}{\partial {\theta }_i}{\mathrm{\Sigma }}^{-1}\frac{\partial \mu }{\partial {\theta }_j}+\frac{1}{{\eta }_{\mathrm{\Sigma }}}\frac{1}{2}tr\left({\mathrm{\Sigma }}^{-1}\frac{\partial \mathrm{\Sigma }}{\partial {\theta }_i}{\mathrm{\Sigma }}^{-1}\frac{\partial \mathrm{\Sigma }}{\partial {\theta }_j}\right).$$

(56)

All of the remainder of this section is simply a rewriting of the work in Section 2 with the twisted Fisher metric instead of the regular Fisher metric. We will use the term “twisted geodesic” instead of “geodesic for the twisted metric”.

This approach seems to be somewhat arbitrary: arguably, the mean and the covariance play a “different role” in the definition of a Gaussian (only the covariance can affect diversity, for example), but we lack a reasonable intrinsic characterization that would make this choice of twisting more natural. This construction can be slightly generalized (see the Appendix).

The IGO flow and the IGO algorithms can be modified to take into account the twisting of the metric; the (η_μ, η_Σ)-twisted IGO flow reads:

$$\frac{\mathrm{d}{\theta }^t}{\mathrm{d}t}=I{({\eta }_{\mu },{\eta }_{\mathrm{\Sigma }})}^{-1}(\theta ){\displaystyle {\int }_X{W}_{{\theta }^t}^f}(x){\nabla }_{\theta }\ln P_{\theta }(x)P_{{\theta }^t}(\mathrm{d}x).$$

(57)

The only difference with (9) is that I⁻¹(θ) has been replaced by I(η_μ, η_Σ)⁻¹(θ).

This leads us to the twisted IGO algorithms.

Definition 12. The (η_μ, η_Σ)-twisted IGO algorithm associated with parametrization θ, sample size N, step size δt and selection scheme w is given by the following update rule:

$${\theta }^{t+\delta t}={\theta }^t+\delta tI{({\eta }_{\mu },{\eta }_{\mathrm{\Sigma }})}^{-1}({\theta }^t){\displaystyle \sum _{i=1}^N{\widehat{w}}_i}\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }.$$

Definition 13. The (η_μ, η_Σ)-twisted geodesic IGO algorithm associated with sample size N, step size δt and selection scheme w is given by the following update rule:

$${\theta }^{t+\delta t}={\exp }_{{\theta }^t}(Y\delta t)$$

(58)

where:

$$Y=I{({\eta }_{\mu },{\eta }_{\mathrm{\Sigma }})}^{-1}({\theta }^t){\displaystyle \sum _{i=1}^N{\widehat{w}}_i\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }}.$$

(59)

By definition, the twisted geodesic IGO algorithm does not depend on the parametrization (but it does depend on η_μ and η_Σ).

There is some redundancy between δt, η_μ and η_Σ: the only values actually appearing in the equations are δtη_μ and δtη_Σ. More formally:

Proposition 9. Let k, d, N ∈ N, η_μ, η_Σ, δt, λ₁, λ₂ ∈ ℝ and w : [0; 1] → ℝ.

The (η_μ, η_Σ)-twisted IGO algorithm with sample size N, step size δt and selection scheme w coincides with the (λ₁η_μ, λ₁η_Σ)-twisted IGO algorithm with sample size N, step size λ₂δt and selection scheme$\frac{1}{{\lambda }_1{\lambda }_2}w$. The same is true for geodesic IGO.

In order to obtain the twisted algorithms, the Fisher metric in IGO has to be replaced by the metric from Definition 11. In practice, the equations found by twisting the metric are exactly the equations without twisting, except that we have “forced” the learning rates η_μ, η_Σ to appear by multiplying the increments of μ and Σ by η_μ and η_Σ.

We can now describe pure rank-μ CMA-ES and xNES with separate learning rates as twisted IGO algorithms:

Proposition 10 (xNES as IGO). The xNES algorithm with sample size N, weights w_i and learning rates η_μ, η_σ = η_B = η_Σ coincides with the$\frac{{\eta }_{\mu }}{\delta t}$, $\frac{{\eta }_{\mathrm{\Sigma }}}{\delta t}$-twisted IGO algorithm with sample size N, weights w_i, step size δt and in which, given the current position (μ_t, A_t), the set of Gaussians is parametrized by:

$$(\delta ,M)\mapsto \mathcal{N}\left({\mu }_t+A_t\delta ,\left(A_t\exp \left(\frac{1}{2}M\right)\right){\left(A_t\exp \left(\frac{1}{2}M\right)\right)}^T\right),$$

with δ ∈ ℝ^m and M ∈ Sym(ℝ^m).

The parameters maintained by the algorithm are (μ, A), and the x_i are sampled from N (μ, AA^T).

Proposition 11 (Pure rank-μ CMA-ES as IGO). The pure rank-μ CMA-ES algorithm with sample size N, weights w_i and learning rates η_μ and η_Σ coincides with the$\left(\frac{{\eta }_{\mu }}{\delta t},\frac{{\eta }_{\mathrm{\Sigma }}}{\delta t}\right)$-twisted IGO algorithm with sample size N, weights w_i, step size δt and the parametrization (μ, Σ).

The proofs of these two statements are an easy rewriting of their non-twisted counterparts: one can return to the non-twisted metric (up to a η_Σ factor) by changing μ to $\frac{\sqrt{\eta \sigma }}{\sqrt{\eta \mu }}\mu$.

We give the equations of the twisted geodesics of ${\mathbb{G}}_d$ in the Appendix.

6.3. Trajectories of Different IGO Steps

As we have seen, two different IGO algorithms (or an IGO algorithm and the GIGO algorithm) coincide at first order in δt when δt → 0. In this section, we study the differences between pure rank-μ CMA-ES, xNES and GIGO by looking at the second order in δt, and in particular, we show that xNES and GIGO do not coincide in the general case.

We view the updates done by one step of the algorithms as paths on the manifold ${\mathbb{G}}_d$, from (μ(t), Σ(t)) to (μ(t + δt), Σ(t + δt)), where δt is the time step of our algorithms, seen as IGO algorithms. More formally:

Definition 14. (1) We call the GIGO update trajectory the application:

$$T_{\mathrm{GIGO}}:(\mu ,\mathrm{\Sigma },v_{\mu }{v}_{\mathrm{\Sigma }})\mapsto (\delta t\mapsto {\exp }_{\mathcal{N}(\mu ,A{A}^T)}(\delta t{\eta }_{\mu }{v}_{\mu },\delta t{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }})).$$

(exp is the exponential of the Riemannian manifold${\mathbb{G}}_d({\eta }_{\mu },{\eta }_{\mathrm{\Sigma }})$)

(2) We call the xNES update trajectory the application:

$$T_{\mathrm{xNES}}:(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})\mapsto (\delta t\mapsto \mathcal{N}(\mu +\delta t{\eta }_{\mu }{v}_{\mu },A\exp [{\eta }_{\mathrm{\Sigma }}\delta t{A}^{-1}{v}_{\mathrm{\Sigma }}{(A^{-1})}^T]A^T)),$$

with AA^T = Σ. The application above does not depend on the choice of a square root A.

(3) We call the CMA-ES update trajectory the application:

$$T_{\mathrm{CMA}}:(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})\mapsto (\delta t\mapsto \mathcal{N}(\mu +\delta t{\eta }_{\mu }{v}_{\mu },A{A}^T+\delta t{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }})).$$

These applications map the set of tangent vectors to${\mathbb{G}}_d$ ($(T{\mathbb{G}}_d)$to the curves in${\mathbb{G}}_d({\eta }_{\mu },{\eta }_{\mathrm{\Sigma }})$.

We will also use the following notation: μ_GIGO := ϕ_μ○T_GIGO, μ_xNES := ϕ_μ○T_xNES, μ_CMA := ϕ_μ○T_CMA, Σ_GIGO := ϕ_Σ ○ T_GIGO, Σ_xNES := ϕ_Σ ○ T_xNES and Σ_CMA := ϕ_Σ ○ T_CMA, where ϕ_μ (resp. ϕ_Σ) extracts the μ-component (resp. the Σ-component) of a curve.

In particular, Im(ϕ_μ) ⊂ ℝ^d and Im(ϕ_Σ) ⊂ P_d, where P_d (the set of real symmetric positive-definitematrices of dimension d) is seen as a subset of ℝ^d2.

For instance, T_GIGO(μ, Σ, v_μ, v_Σ)(δt) gives the position (mean and covariance matrix) of the GIGO algorithm after a step of size δt, while μ_GIGO and Σ_GIGO give, respectively, the mean component and the covariance component of this position.

This formulation ensures that the trajectories we are comparing had the same initial position and the same initial speed, which is the case provided the sampled points (the values directly sampled from $\mathcal{N}(\mu ,\mathrm{\Sigma })$, not from $\mathcal{N}(0,I)$ and transformed) are the same.

Different IGO algorithms coincide at first order in δt. The following proposition gives the second order expansion of the trajectories of the algorithms.

Proposition 12 (Second derivatives of the trajectories). We have:

$$\begin{array}{c}{\mu }_{\mathrm{GIGO}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)={\eta }_{\mu }{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }}{\displaystyle {\sum }_0^{-1}{v}_{\mu },}\\ {\mu }_{\mathrm{xNES}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)={\mu }_{\mathrm{CMA}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)=0,\\ {\mathrm{\Sigma }}_{\mathrm{GIGO}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)={\eta }_{\mathrm{\Sigma }}^2{v}_{\mathrm{\Sigma }}{\mathrm{\Sigma }}^{-1}{v}_{\mathrm{\Sigma }}-{\eta }_{\mu }{\eta }_{\mathrm{\Sigma }}{v}_{\mu }{v}_{\mu }^T,\\ {\mathrm{\Sigma }}_{\mathrm{xNES}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)={\mu }_{\mathrm{\Sigma }}^2{v}_{\mathrm{\Sigma }}{\mathrm{\Sigma }}^{-1}{v}_{\mathrm{\Sigma }},\\ {\mathrm{\Sigma }}_{\mathrm{CMA}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)=0.\end{array}$$

Proof. We can immediately see that the second derivatives of μ_xNES, μ_CMA and Σ_CMA are zero. Next, we have:

$$\begin{array}{l}{\mathrm{\Sigma }}_{\mathrm{xNES}}(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})(t)=A\exp [t{A}^{-1}{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }}{(A^{-1})}^T]A^T\\ =A{A}^T+t{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }}+\frac{t^2}{2}{\eta }_{\mathrm{\Sigma }}^2{v}_{\mathrm{\Sigma }}{(A^{-1})}^T{A}^{-1}{v}_{\mathrm{\Sigma }}+o(t^2)\\ =\mathrm{\Sigma }+t{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }}+\frac{t^2}{2}{\eta }_{\mathrm{\Sigma }}^2{v}_{\mathrm{\Sigma }}{\mathrm{\Sigma }}^{-1}{v}_{\mathrm{\Sigma }}+o(t^2).\end{array}$$

The expression of Σ_xNES(μ, Σ, v_μ, v_Σ)^″(0) follows.

Now, for GIGO, let us consider the geodesic starting at (μ₀, Σ₀) with initial speed (η_μv_μ, η_Σv_Σ). By writing J_μ(0) = J_μ(t), we find $\dot{\mu }(t)=\mathrm{\Sigma }(t){\mathrm{\Sigma }}_0^{-1}{\dot{\mu }}_0$. We then easily have $\ddot{\mu }(0)={\dot{\mathrm{\Sigma }}}_0{\mathrm{\Sigma }}_0^{-1}{\dot{\mu }}_0$ In other words:

$${\mu }_{\mathrm{GIGO}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)={\eta }_{\mu }{\eta }_{\mathrm{\Sigma }}{v}_{\mathrm{\Sigma }}{\displaystyle {\sum }_0^{-1}{v}_{\mu }.}$$

Finally, by using Theorem 4 and differentiating, we find:

$$\begin{array}{c}\ddot{\mathrm{\Sigma }}={\eta }_{\mathrm{\Sigma }}\dot{\mathrm{\Sigma }}(J_{\mathrm{\Sigma }}-J_{\mu }{\mu }^T)-{\eta }_{\mathrm{\Sigma }}\mathrm{\Sigma }{J}_{\mu }{\dot{\mu }}^T,\\ {\ddot{\mathrm{\Sigma }}}_0={\eta }_{\mathrm{\Sigma }}{\dot{\mathrm{\Sigma }}}_0\frac{1}{{\eta }_{\mathrm{\Sigma }}}{\mathrm{\Sigma }}_0^{-1}{\dot{\mathrm{\Sigma }}}_0-\frac{{\eta }_{\mathrm{\Sigma }}}{{\eta }_{{}^{{}^{\mu }}}}{\dot{\mu }}_0{\dot{\mu }}_0^T={\eta }_{\mathrm{\Sigma }}^2{v}_{\mathrm{\Sigma }}{\mathrm{\Sigma }}_0^{-1}{v}_{\mathrm{\Sigma }}-{\eta }_{\mathrm{\Sigma }}{\eta }_{\mu }{v}_{\mu }{v}_{\mu }^T.\end{array}$$

In order to interpret these results, we will look at what happens in dimension one. In higher dimensions, we can suppose that the algorithms exhibit a similar behavior, but an exact interpretation is more difficult for GIGO in ${\mathbb{G}}_d$.

• In [19], it has been noted that xNES converges to quadratic minima slower than CMA-ES and that it is less subject to premature convergence. That fact can be explained by observing that the mean update is exactly the same for CMA-ES and xNES, whereas xNES tends to have a higher variance (Proposition 12 shows this at order two, and it is easy to see that in dimension one, for any μ, Σ, v_μ, v_Σ, we have Σ_xNES(μ, Σ, v_μ, v_Σ) > Σ_CMA(μ, Σ, v_μ, v_Σ)).
• At order two, GIGO moves the mean faster than xNES and CMA-ES if the standard deviation is increasing and more slowly if it is decreasing. This seems to be a reasonable behavior (if the covariance is decreasing, then the algorithm is presumably close to a minimum, and it should not leave the area too quickly). This remark holds only for isolated steps, because we do not take into account the evolution of the variance.
• The geodesics of ${\mathbb{G}}_1$ are half-circles (see Figure 2 below; we recall that ${\mathbb{G}}_1$ is the Poincaré half-plane). Consequently, if the mean is supposed to move (which always happens), then σ → 0 when δt → ∞. For example, a step whose initial speed has no component on the standard deviation will always decrease it. See also Proposition 15, about the optimization of a linear function.
• For the same reason, for a given initial speed, the update of μ always stays bounded as a function of δt: it is not possible to make one step of the GIGO algorithm go further than a fixed point by increasing δt. Still, the geodesic followed by GIGO changes at each step, so the mean of the overall algorithm is not bounded.

We now show that xNES follows the geodesics of ${\mathbb{G}}_d$ if the mean is fixed, but that xNES and GIGO do not coincide otherwise.

Proposition 13 (xNES is not GIGO in the general case). Let μ, v_μ ∈ ℝ^d, A ∈ GL_d, v_Σ ∈ M_d.

Then, the GIGO and xNES updates starting at$\mathcal{N}(\mu ,\mathrm{\Sigma })$ with initial speeds v_μ and v_Σ follow the same trajectory if and only if the mean remains constant. In other words:

$$T_{\mathrm{GIGO}}(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})=T_{\mathrm{xNES}}(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})ifandonlyif{v}_{\mu }=0.$$

Proof. If v_μ = 0, then we can compute the GIGO update by using Theorem 4: since J_μ = 0, $\dot{\mu }=0$, and μ remains constant. Now, we have $J_{\mathrm{\Sigma }}={\mathrm{\Sigma }}^{-1}\dot{\mathrm{\Sigma }}$; this is enough information to compute the update. Since this quantity is also preserved by the xNES algorithm (see, for example, the proof of Proposition 14), the two updates coincide.

Figure 2. One step of the geodesic IGO (GIGO) update.

If v_μ ≠ 0, then ${\mathrm{\Sigma }}_{\mathrm{xNES}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)-{\mathrm{\Sigma }}_{\mathrm{GIGO}}{(\mu ,\mathrm{\Sigma },v_{\mu },v_{\mathrm{\Sigma }})}^{″}(0)={\eta }_{\mu }{\eta }_{\mathrm{\Sigma }}{v}_{\mu }{v}_{\mu }^T\ne 0$ and, in particular, T_GIGO(μ, Σ, v_μ, v_Σ) ≠ T_xNES(μ, Σ, v_μ, v_Σ).

6.4. Blockwise GIGO

Although xNES is not GIGO, it is possible to define a family of algorithms extending GIGO and including xNES, by decomposing our family of probability distributions as a product and by following the restricted geodesics simultaneously.

Definition 15 (Splitting). Let Θ be a Riemannian manifold. A splitting of Θ is n manifolds Θ₁, …, Θ_n and a diffeomorphism Θ ≅ Θ₁ × … × Θ_n. If for all x ∈ Θ, for all 1 ≤ i < j ≤ n, we also have T_i,xM ⊥ T_j,xM as subspaces of T_xM (see Notation 2), then the splitting is said to be compatible with the Riemannian structure. If the Riemannian manifold is not ambiguous, we will simply write a “compatible splitting”.

We now give some notation, and we define the blockwise GIGO update:

Notation 2. Let Θ be a Riemannian manifold, Θ₁, …, Θ_n a splitting of Θ, θ = (θ₁, …, θ_n) ∈ Θ, Y ∈ T_θΘ and 1 ≤ i ≤ n.

• We denote by Θ_θ,i the Riemannian manifold

  $$\left\{{\theta }_1\right\}\times \dots \times \left\{{\theta }_{i-1}\right\}\times {\mathrm{\Theta }}_i\times \left\{{\theta }_{i+1}\right\}\times \dots \times \left\{{\theta }_n\right\},$$

  with the metric induced from Θ. There is a canonical isomorphism of vector spaces$T_{\theta }\mathrm{\Theta }={\oplus }_{i=1}^{n}T{\mathrm{\Theta }}_{\theta ,i}$. Moreover, if the splitting is compatible, it is an isomorphism of Euclidean spaces.
• We denote by Φ_θ,i the exponential at θ of the manifold Θ_θ,i.

Definition 16 (Blockwise GIGO update). Let Θ₁, …, Θ_n be a compatible splitting. The blockwise GIGO algorithm in Θ with splitting Θ₁, …, Θ_n associated with sample size N, step sizes δt₁, …, δt_n and selection scheme w is given by the following update rule:

$$\theta \leftarrow ({\theta }_1^{t+\delta t_1},\dots ,{\theta }_n^{t+\delta t_n})$$

(60)

where:

$$Y=I^{-1}({\theta }^t){\displaystyle \sum _{i=1}^N{\widehat{w}}_i\frac{\partial \ln P_{\theta }(x_i)}{\partial \theta }},$$

(61)

$${\theta }_k^{t+\delta t_k}={\mathrm{\Phi }}_{{\theta }^t,k}(\delta t_k{Y}_k),$$

(62)

with Y_k the TΘ_θ,k-component of Y. This update only depends on the splitting (and not on the parametrization inside each Θ_k).

The compatibility condition ensures that the natural gradient of $W_{{\theta }^t}^f$ (defined in Section 2.2) in the whole manifold Θ really is the sum of the gradients of this same function in the submanifolds Θ_k. A practical consequence is that the Y_k in Equation (62) can be computed simply by taking the natural gradient in Θ_k:

$$Y_k=I_k^{-1}({\theta }_i^t){\displaystyle \sum _{i=1}^N\widehat{w}}\frac{\partial \ln P_{\theta }(x_i)}{\partial {\theta }_k},$$

(63)

where I_k is the metric of Θ_k.

Since blockwise GIGO only depends on the splitting (and the tunable parameters: sample size, step sizes and selection scheme), it can be thought of as almost parametrization-invariant.

Notice that blockwise GIGO updates and twisted GIGO updates are two different things: firstly, blockwise GIGO can be defined on any manifold with a compatible splitting, whereas twisted GIGO (and twisted IGO) are only defined for Gaussians. However, even in ${\mathbb{G}}_d({\eta }_{\mu },{\eta }_{\sum })$, with the splitting (μ, Σ), these two algorithms are different: for instance, if η_μ = η_Σ and δt = 1, then the twisted GIGO is the regular GIGO algorithm, whereas blockwise GIGO is not (actually, we will prove that it is the xNES algorithm). The only thing blockwise GIGO and twisted GIGO have in common is that they are compatible with the (η_μ, η_Σ)-twisted IGO flow Equation (57): a parameter θ^t following these updates with δt → 0 and N → ∞ is a solution of Equation (57).

We now have a new description of the xNES algorithm:

Proposition 14 (xNES is a Blockwise GIGO algorithm). The Blockwise GIGO algorithm in${\mathbb{G}}_d$ with splitting$\mathrm{\Phi }:\mathcal{N}\left(\mu ,\mathrm{\Sigma }\right)\mapsto \left(\mu ,\mathrm{\Sigma }\right)$, sample size N, step sizes δt_μ, δt_Σ and selection scheme w coincides with the xNES algorithm with sample size N, weights w_i and learning rates η_μ = δt_μ, η_σ = η_B = δt_Σ.

Proof. Firstly, notice that the splitting (μ, Σ) is compatible, by Proposition 1.

Now, let us compute the Blockwise GIGO update: we have ${\mathbb{G}}_d\cong {ℝ}^d\times P_d$, where P_d is the space of real positive-definite matrices of dimension d. We have ${\mathrm{\Theta }}_{{\theta }^t,1}={ℝ}^d\times (\{{\sum }^t\})\to {\mathbb{G}}_d,{\mathrm{\Theta }}_{{\theta }^t,2}=(\{{\mu }^t\}\times P_d)\to {\mathbb{G}}_d$. The induced metric on ${\mathrm{\Theta }}_{{\theta }^t}{}_{{,}_1}$ is the Euclidean metric, so we have:

$$\mu \leftarrow {\mu }^t+\delta t_1{Y}_{\mu }.$$

Since we have already shown (using the notation in Definition 9) that Y_μ = AG_μ (in the proof of Proposition 6), we find:

$$\mu \leftarrow {\mu }^t+\delta t_{1}A{\mathbb{G}}_{\mu }.$$

On Θ_θt,₂, we have the following Lagrangian for the geodesics:

$$ℒ(\sum ,\dot{\sum })=\frac{1}{2}tr(\dot{\sum }{\sum }^{-1}\dot{\sum }{\sum }^{-1}).$$

By applying Noether’s theorem, we find that

$$J_{\mathrm{\Sigma }}={\mathrm{\Sigma }}^{-}{}^1\dot{\mathrm{\Sigma }}$$

is invariant along the geodesics of ${\mathrm{\Theta }}_{{\theta }^t,2}$, so they are defined by the equation $\dot{\mathrm{\Sigma }}=\mathrm{\Sigma }{J}_{\mathrm{\Sigma }}=\mathrm{\Sigma }{\mathrm{\Sigma }}_0^{-1}{\dot{\mathrm{\Sigma }}}_0$ (and therefore, any update preserving the invariant J_Σ will satisfy this first-order differential equation and follow the geodesics of ${\mathrm{\Theta }}_{{\theta }^t,2})$. The xNES update for the covariance matrix is given by A(t) = A₀ exp(tG_M/2). Therefore, we have $\mathrm{\Sigma }\left(t\right)=A_0\exp \left(t{\mathbb{G}}_M\right)A_0^T,$${\mathrm{\Sigma }}^{-1}(t)={(A_0^{-1})}^T\exp \left(-t{\mathbb{G}}_M\right)A_0^{-1}$, $\dot{\sum }(t)=A_0\exp (t{\mathbb{G}}_M){\mathbb{G}}_M{A}_0^T$ and, finally, ${\mathrm{\Sigma }}^{-}{}^1(t)\dot{\mathrm{\Sigma }}(t)={(A_0^{-1})}^T{\mathbb{G}}_M{A}_0^T={\sum }_0^{-1}{\dot{\sum }}_0$. Therefore, xNES preserves J_Σ, and therefore, xNES follows the geodesics of ${\mathrm{\Theta }}_{{\theta }^t,2}$ (notice that we had already proven this in Proposition 13, since we are looking at the geodesics of ${\mathbb{G}}_d$ with a fixed mean).

Although blockwise GIGO is somewhat “less natural” than GIGO, it can be easier to compute for some splittings (as we have just seen), and in the case of the Gaussian distributions, the mean-covariance splitting seems reasonable.

7. Numerical Experiments

We conclude this article with some numerical experiments to compare the behavior of GIGO, xNES and pure rank-μ CMA-ES (we give the pseudocodes for these algorithms in the Appendix). We made two series of tests. The first one is a performance test, using classical benchmark functions and the settings from [19]. The goal of the second series of tests is to illustrate the computations in Section 6.3 by plotting the trajectories (standard deviation versus mean) of these three algorithms in dimension one.

The source code is available at [22].

7.1. Benchmarking

For the first series of experiments, presented in Figure 3, we used the following parameters, taken from [19] (we recall that xNES and pure rank-μ CMA-ES are seen as IGO algorithms):

• Varying dimension.
• Sample size: $\lfloor 4+3\log (d)\rfloor .$
• Weights: $w_i=\frac{max\left(0,\log \left(\frac{n}{2}+1\right)-\log (i)\right)}{{\sum }_{j=1}^{N}max\left(0,\log \left(\frac{n}{2}+1\right)-\log (j)\right)}-\frac{1}{N}$.
• IGO step size and learning rates: δt = 1, η_μ = 1, ${\eta }_{\sum }=\frac{3}{5}\frac{3+\log (d)}{d\sqrt{d}}$..
• Initial position: ${\theta }^0=\mathcal{N}\left(x_0,I\right)$, where x₀ is a random point of the circle with center zero, and radius 10.
• Euler method for GIGO: Number of steps: 100. We used the GIGO-A variant of the algorithm. No significant difference was noticed with GIGO-Σ or with the exact GIGO algorithm. The only advantage of having an explicit solution of the geodesic equations is that the update is quicker to compute.
• We chose not to use the exact expression of the geodesics for this benchmarking to show that having to use the Euler method is fine. However, we did run the tests, and the results are basically the same as GIGO-A.

We plot the median number of runs to achieve target fitness (10⁻⁸). Each algorithm has been tested in dimension 2, 4, 8, 16, 32 and 64: a missing point means that all runs converged prematurely.

7.1.1. Failed Runs

In Figure 3, a point is plotted even if only one run was successful. Below is the list of the settings for which at least one run converged prematurely.

• Only one run reached the optimum for the cigar-tablet function with CMA-ES in dimension eight.
• Seven runs (out of 24) reached the optimum for the Rosenbrock function with CMA-ES in dimension 16.
• About half of the runs reached the optimum for the sphere function with CMA-ES in dimension four.

For the following settings, all runs converged prematurely.

• GIGO did not find the optimum of the Rosenbrock function in any dimension.
• CMA-ES did not find the optimum of the Rosenbrock function in dimension 2, 4, 32 and 64.
• All of the runs converged prematurely for the cigar-tablet function in dimension two with CMA-ES, for the sphere function in dimension two for all algorithms and for the Rosenbrock function in dimension two and four for all algorithms.

7.1.2. Discussion

As the last item in Section 7.1.1 shows, all of the algorithms converge prematurely in a low dimension, probably because the covariance learning rate has been set too high (or because the sample size is too small). This is different from the results in [19].

This remark aside, as noted in [19], the xNES algorithm shows more robustness than CMA-ES and GIGO: it is the only algorithm able to find the minimum of the Rosenbrock function in high dimensions. However, its convergence is consistently slower.

In terms of performance, when both of them work, pure rank-μ CMA-ES (or equivalently, IGO in the parametrization (μ, Σ)) and GIGO are extremely close (GIGO is usually a bit better). An advantage of GIGO is that it is theoretically defined for any δt, η_Σ, whereas the covariance matrix maintained by CMA-ES (not only pure rank-μ CMA-ES) can stop being positive definite if η_Σδt > 1. However, in that case, the GIGO algorithm is prone to premature convergence (remember Figure 2 and see Proposition 15 below), and in practice, the learning rates are much smaller.

7.2. Plotting Trajectories in ${\mathbb{G}}_1$

We want the second series of experiments to illustrate the remarks about the trajectories of the algorithms in Section 6.3, so we decided to take a large sample size to limit randomness, and we chose a fixed starting point for the same reason. We use the weights below because of the property of quantile improvement proven in [23]: the 1/4-quantile will improve at each step. The parameters we used were the following:

• Sample size: λ = 5, 000
• Dimension one only.
• Weights: w = 41_q⩽1/4 (w_i = 4.1_i⩽1,250)
• IGO step size and learning rates: η_μ = 1, ${\eta }_{\sum }=\frac{3}{5}\frac{3+\log (d)}{d\sqrt{d}}=1.8$, varying δt.
• Initial position: ${\theta }^0=\mathcal{N}\left(10,1\right)$
• Dots are placed at t = 0, 1, 2 … (except for the graph δt = 1.5, for which there is a dot for each step).

Figures 4–8 show the optimization of x ↦ x², and Figures 9–11 show the optimization of x ↦ −x.

Figures 7, 8 and 11 show that when δt ≥ 1, GIGO reduces the covariance, even at the first step. More generally, when using the GIGO algorithm in ${\tilde{\mathbb{G}}}_d$ for the optimization of a linear function, there exists a critical step size δt_cr (depending on the learning rates η_μ, η_σ and on the weights w_i), above which, GIGO will converge, and we can compute its value when the weights are of the form $1_q{}_{\le q_0}$ (for q₀ ≥ 0.5, the discussion is not relevant, because in that case, even the IGO flow converges prematurely. Compare with the critical δt of the smoothed cross entropy method and IGO-ML in [1]).

Proposition 15. Let d ∈ ℕ, k, η_μ, ${\eta }_{\sigma }\in {ℝ}_{+}^{*}$; let$w=k{.1}_{q\le q_0}$ and let

$$\begin{array}{r}\hfill \begin{array}{rrrrr}\hfill g& \hfill :& \hfill {ℝ}^d& \hfill \to & \hfill ℝ\end{array}\\ \hfill \begin{array}{rrr}\hfill x& \hfill \mapsto & \hfill -x_1\end{array}\end{array}.$$

Let μ_n be the first coordinate of the mean, and let${\sigma }_n^2$ be the variance (at step n) maintained by the (η_μ, η_σ)-twisted geodesic IGO algorithm in${\tilde{\mathbb{G}}}_d$ associated with selection scheme w, sample size ∞ and step size δt, when optimizing g (“sample size ∞” meaning the limit of the update when the sample size tends to infinity, which is deterministic [1]).

There exists δt_cr, such that:

• if δt > δt_cr, (σ_n) converges to zero with exponential speed and (μ_n) converges.
• if δt = δt_cr, (σ_n) remains constant and (μ_n) tends to ∞ with linear speed.
• if 0 < δt < δt_cr, both (σ_n) and μ_n tend to ∞ with exponential speed.

The proof and the expression of δt_cr can be found in the Appendix.

In the case corresponding to k = 4, n = 1, q₀ = 1/4, η_μ = 1, η_σ = 1.8, we find:

$$\delta t_{\mathrm{cr}}\approx \mathrm{0.84.}$$

(64)

8. Conclusions

We introduced the geodesic IGO algorithm, and we showed that in the case of Gaussian distributions, Noether’s theorem directly gives a first order equation satisfied by the geodesics. In terms of performance, the GIGO algorithm is similar to pure rank-μ CMA-ES, which is rather encouraging: it would be interesting to test GIGO on real problems. Moreover, GIGO is a reasonable and totally parametrization-invariant algorithm (provided we can compute the solution of the equations of the geodesics), and as such, it should be studied for other families of probability distributions, like Bernoulli distributions (although in that case, the Riemannian exponential is not defined if the step size is too large, because the length of the geodesics is finite). Noether’s theorem could be a crucial tool for this.

We also showed that xNES and GIGO are not the same algorithm, and we defined blockwise GIGO, a simple extension of the GIGO algorithm, showing that xNES has a special status, as it admits a definition that is “almost” parametrization-invariant.