2. Some examples🔗

In this page, we present two examples of algorithms encompassed by our framework: the Pure Random Search (PRS) and the LIPO algorithms. The formalization of these algorithms in our framework relies on the definition of the initial probability measure and the Markov kernels that define how to sample the next element based on the previous ones.

2.1. Pure Random Search🔗

The Pure Random Search (PRS) algorithm is a simple stochastic iterative global optimization algorithm that samples uniformly from the input space at each iteration using a fixed probability measure μ. It can be represented in our framework as follows:

def

PRS.{u_1, u_2} {α : Type u_1} {β : Type u_2} [MeasurableSpace α]
  [MeasurableSpace β] (μ : MeasureTheory.Measure α)
  [MeasureTheory.IsProbabilityMeasure μ] : Algorithm α β
PRS.{u_1, u_2} {α : Type u_1}
  {β : Type u_2} [MeasurableSpace α]
  [MeasurableSpace β]
  (μ : MeasureTheory.Measure α)
  [MeasureTheory.IsProbabilityMeasure μ] :
  Algorithm α β

The Pure Random Search (PRS) algorithm for global optimization. This baseline algorithm samples uniformly from the input space at each iteration using a fixed probability measure μ.

noncomputable def PRS : Algorithm α β where
  ν := μ
  kernel_iter _ := Kernel.const _ μ

2.2. Decision-based methods🔗

We provide an interface for decision-based optimization algorithms. Such algorithm are described by a sequence of decision rules to decide wether to accept or reject a candidate solution at each iteration, based on the observed data. Decision rules are functions with the signature D_n : \Omega^{n + 1} \times \beta^{n + 1} \to [0, 1], where \Omega is the search space and \beta is the space of function values (most likely \mathbb{R}). More formally, given a function to optimize V : \Omega \to \beta at each iteration n > 0, the algorithm samples a candidate x from a distribution \mu on \Omega and accepts it with probability D_n(C_n, V(C_n), x).

We define Decision as a special case of the Algorithm structure, for which the Markov kernel at iteration n is defined through the decision rules (D_n)_{n \in \mathbb{N}}. Measurability assumptions are required on the decision rules for the kernels to be well-defined and we only support \{0, 1\}-valued decision rules at the moment.

def

Decision.{u_1, u_2} {α : Type u_1} {β : Type u_2} [MeasurableSpace α]
  [MeasurableSpace β] (μ : MeasureTheory.Measure α)
  [MeasureTheory.IsProbabilityMeasure μ]
  {decision : (n : ℕ) → prod_iter_image α β n → Set α}
  (measurableSet_decision_prod :
    ∀ (n : ℕ), MeasurableSet {p | p.2 ∈ decision n p.1})
  (h :
    ∀ (n : ℕ) (data : prod_iter_image α β n), μ (decision n data) ≠ 0) :
  Algorithm α β
Decision.{u_1, u_2} {α : Type u_1}
  {β : Type u_2} [MeasurableSpace α]
  [MeasurableSpace β]
  (μ : MeasureTheory.Measure α)
  [MeasureTheory.IsProbabilityMeasure μ]
  {decision :
    (n : ℕ) →
      prod_iter_image α β n → Set α}
  (measurableSet_decision_prod :
    ∀ (n : ℕ),
      MeasurableSet
        {p | p.2 ∈ decision n p.1})
  (h :
    ∀ (n : ℕ)
      (data : prod_iter_image α β n),
      μ (decision n data) ≠ 0) :
  Algorithm α β

The interface for decision-based optimization algorithms.

2.3. LIPO🔗

The LIPO algorithm has been introduced in (Malherbe and Vayatis, 2017)Cédric Malherbe and Nicolas Vayatis, 2017. “Global optimization of Lipschitz functions”. In International conference on machine learning. and is a decision-based global optimization algorithm made for Lipschitz functions. It uses the Lipschitz constant to adaptively construct an upper bound that guides the sampling of the search space. More formally, at each iteration, LIPO samples from the set of potential maximizers defined as \left \{ x \in \alpha \; \middle| \; \max_{1 \le i \le n} f(X_i) \le \min_{1 \le i \le n} f(X_i) + \kappa \|X_i - x \|_2 \right\}, where α is a pseudo-metric, measurable, Borel, and second-countable space and \kappa is the Lipschitz constant of the function f : \alpha \to \mathbb{R}.

It can be represented in our framework as a special case of the Decision structure, where the decision rules return 1 if the candidate solution is in the set of potential maximizers and 0 otherwise:

noncomputable def LIPO : Algorithm α ℝ :=
  Decision μ (fun _ ↦ measurableSet_potential_max_prod κ) h

Note that LIPO requires the set of potential maximizers to have non-zero measure at each iteration, ensuring that the algorithm can sample from it. This is a non-trivial assumption that depends on the choice of the initial probability measure μ, the function to optimize, and the \sigma-algebra on the search space.

def

LIPO.{u_1} {α : Type u_1} [PseudoMetricSpace α] [MeasurableSpace α]
  [BorelSpace α] [SecondCountableTopology α]
  (μ : MeasureTheory.Measure α) [MeasureTheory.IsProbabilityMeasure μ]
  (κ : NNReal)
  (h :
    ∀ (n : ℕ) (data : prod_iter_image α ℝ n),
      μ (LIPO.potential_max κ data) ≠ 0) :
  Algorithm α ℝ
LIPO.{u_1} {α : Type u_1}
  [PseudoMetricSpace α]
  [MeasurableSpace α] [BorelSpace α]
  [SecondCountableTopology α]
  (μ : MeasureTheory.Measure α)
  [MeasureTheory.IsProbabilityMeasure μ]
  (κ : NNReal)
  (h :
    ∀ (n : ℕ)
      (data : prod_iter_image α ℝ n),
      μ (LIPO.potential_max κ data) ≠ 0) :
  Algorithm α ℝ

The LIPO (LIPschitz Optimization) algorithm for global optimization. This algorithm optimizes an unknown function assuming only that it has a finite Lipschitz constant κ. It starts with an arbitrary probability measure μ as initial distribution and iteratively samples from the set of potential maximizers, ensuring consistency and convergence to the global optimum (Malherbe et al., 2017).

2.4. RankOpt🔗

The RankOpt algorithm has been introduced in (Malherbe and Vayatis, 2017)Cédric Malherbe and Nicolas Vayatis, 2017. “A Ranking Approach to Global Optimization”. In International conference on machine learning. and is a sophisticated decision-based global optimization algorithm. It is based on the notion of a ranking rule, which is a function induced by another function f and is defined as r_f (x, y) = \begin{cases} 1 & \text{if } f(x) > f(y) \\ 0 & \text{if } f(x) = f(y) \\ -1 & \text{if } f(x) < f(y) \end{cases}. Ranking rules define equivalence classes of functions that share the same induced ranking. To use RankOpt, one needs to define a countable set of ranking rules \mathcal{R} (e.g. a subset of the ranking rules induced by continuous functions). The algorithm samples points x such that, for any ranking rule r \in \mathcal{R} that is consistent with the observed data, 0 \le r(x, \arg \max_{1 \le i \le n} f(X_i)). The set of such points is the set of potential maximizers of the algorithm and is defined as \left \{x \in \alpha \; \middle| \; \exists r \in \mathcal{R}, \; \mathcal{L}_n(r) = 0 \land 0 \le r(x, \arg \max_{1 \le i \le n} f(X_i))\right\}, where α is a measurable space and \mathcal{L}_n(r) is the ranking loss of r on the observed data, defined as \mathcal{L}_n(r) \triangleq \frac{2}{n (n + 1)} \sum_{1 \le i \le j \le n} \mathbb{I}\left[r(X_i, X_j) \neq r_f(X_i, X_j)\right]. Note that RankOpt does not require to know r_f explicitly as it is evaluated only on observed data points.

noncomputable def RankOpt : Algorithm α β :=
  Decision μ (fun _ ↦ measurableSet_potential_max_prod h𝓡) h

Note that RankOpt requires the set of potential maximizers to have non-zero measure at each iteration, ensuring that the algorithm can sample from it. This is a non-trivial assumption that depends on the choice of the initial probability measure μ, the function to optimize, and the \sigma-algebra on the search space.

def

RankOpt.{u_1, u_2} {α : Type u_1} {β : Type u_2} [MeasurableSpace α]
  (μ : MeasureTheory.Measure α) [MeasureTheory.IsProbabilityMeasure μ]
  [TopologicalSpace β] [MeasurableSpace β] [BorelSpace β]
  [LinearOrder β] [SecondCountableTopology β] [OpensMeasurableSpace β]
  [OrderClosedTopology β] {𝓡 : Set (RankRule α)} (h𝓡 : 𝓡.Countable)
  (h :
    ∀ (n : ℕ) (data : prod_iter_image α β n),
      μ (RankOpt.potential_max data 𝓡) ≠ 0) :
  Algorithm α β
RankOpt.{u_1, u_2} {α : Type u_1}
  {β : Type u_2} [MeasurableSpace α]
  (μ : MeasureTheory.Measure α)
  [MeasureTheory.IsProbabilityMeasure μ]
  [TopologicalSpace β] [MeasurableSpace β]
  [BorelSpace β] [LinearOrder β]
  [SecondCountableTopology β]
  [OpensMeasurableSpace β]
  [OrderClosedTopology β]
  {𝓡 : Set (RankRule α)}
  (h𝓡 : 𝓡.Countable)
  (h :
    ∀ (n : ℕ)
      (data : prod_iter_image α β n),
      μ (RankOpt.potential_max data 𝓡) ≠
        0) :
  Algorithm α β

The RankOpt algorithm for global optimization. This algorithm uses a ranking approach to optimize an unknown function. It maintains a hypothesis class 𝓡 of ranking rules. It starts with an arbitrary probability measure μ as initial distribution and samples from the set of points that could be optimal according to ranking rules consistent with the observed data (Malherbe et al., 2017).

The term RankRule is defined as the subtype of \{-1, 0, 1\}-valued functions that are jointly measurable:

def RankRule (α : Type*) [MeasurableSpace α] :=
  {f : α → α → ({-1, 0, 1} : Set ℝ) // Measurable <| Function.uncurry f}

2.5. CMA-ES🔗

A general implementation of the CMA-ES algorithm in any dimension. As CMA-ES samples \lambda points at each iteration, the input space of the algorithm is \mathbb{R}^{d \times \lambda}, which represents a sequence of \lambda points in \mathbb{R}^{d}. The initial measure is the product of \lambda standard multivariate Gaussian measures on \mathbb{R}^{d}, and the kernel is defined as a product of \lambda multivariate Gaussian measures, where the mean and covariance matrix are given by measurable functions of the past evaluations. These functions can be anything as long as they are measurable w.r.t. the history of the algorithm, thus allowing for any CMA-ES variant to be implemented in this framework.

def

CMA_ES.{u_1} (d lam : ℕ) {β : Type u_1} [MeasurableSpace β]
  {mean :
    (n : ℕ) → prod_iter_image (ℝ_ d lam) β n → EuclideanSpace ℝ (Fin d)}
  (hmean : ∀ (n : ℕ), Measurable (mean n))
  {covar :
    (n : ℕ) → prod_iter_image (ℝ_ d lam) β n → Matrix (Fin d) (Fin d) ℝ}
  (hcovar : ∀ (n : ℕ), Measurable (covar n))
  (m : EuclideanSpace ℝ (Fin d)) (S : Matrix (Fin d) (Fin d) ℝ) :
  Algorithm (ℝ_ d lam) β
CMA_ES.{u_1} (d lam : ℕ) {β : Type u_1}
  [MeasurableSpace β]
  {mean :
    (n : ℕ) →
      prod_iter_image (ℝ_ d lam) β n →
        EuclideanSpace ℝ (Fin d)}
  (hmean : ∀ (n : ℕ), Measurable (mean n))
  {covar :
    (n : ℕ) →
      prod_iter_image (ℝ_ d lam) β n →
        Matrix (Fin d) (Fin d) ℝ}
  (hcovar :
    ∀ (n : ℕ), Measurable (covar n))
  (m : EuclideanSpace ℝ (Fin d))
  (S : Matrix (Fin d) (Fin d) ℝ) :
  Algorithm (ℝ_ d lam) β

The Covariance Matrix Adaptation - Evolution Strategy (CMA-ES) algorithm for global optimization, given the mean and covariance update rules as measurable functions of the history.

noncomputable def CMA_ES : Algorithm (ℝ_ d lam) β where
  ν := Measure.pi (fun _ ↦ multivariateGaussian m S)
  kernel_iter := CMAKernel d lam hmean hcovar
  markov_kernel n := ⟨fun a => byd:ℕlam:ℕβ:Type u_1inst✝:MeasurableSpace βmean:(n : ℕ) → prod_iter_image (ℝ_ d lam) β n → EuclideanSpace ℝ (Fin d)hmean:∀ (n : ℕ), Measurable (mean n)covar:(n : ℕ) → prod_iter_image (ℝ_ d lam) β n → Matrix (Fin d) (Fin d) ℝhcovar:∀ (n : ℕ), Measurable (covar n)m:EuclideanSpace ℝ (Fin d)S:Matrix (Fin d) (Fin d) ℝn:ℕa:prod_iter_image (ℝ_ d lam) β n⊢ IsProbabilityMeasure ((CMAKernel d lam hmean hcovar n) a) simp [CMAKernel]d:ℕlam:ℕβ:Type u_1inst✝:MeasurableSpace βmean:(n : ℕ) → prod_iter_image (ℝ_ d lam) β n → EuclideanSpace ℝ (Fin d)hmean:∀ (n : ℕ), Measurable (mean n)covar:(n : ℕ) → prod_iter_image (ℝ_ d lam) β n → Matrix (Fin d) (Fin d) ℝhcovar:∀ (n : ℕ), Measurable (covar n)m:EuclideanSpace ℝ (Fin d)S:Matrix (Fin d) (Fin d) ℝn:ℕa:prod_iter_image (ℝ_ d lam) β n⊢ IsProbabilityMeasure (Measure.pi fun x => multivariateGaussian (mean n a) (covar n a)); infer_instanceAll goals completed! 🐙⟩

2.5.1. The original CMA-ES🔗

The historical instantiation of this scheme (Hansen and Ostermeier, 1996)Nikolaus Hansen and Andreas Ostermeier, 1996. “Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation”. In Proceedings of IEEE international conference on evolutionary computation. ranks the \lambda points of each generation according to their evaluations and recombines the \mu best ones. It adapts a state made of the mean m, the step size \sigma, the covariance matrix C and two evolution paths p_c and p_\sigma, the points of a generation being sampled i.i.d. according to \mathcal{N}(m, \sigma^2 C).

2.5.1.1. Ranking a generation🔗

As LeanGO maximizes objective functions, a point of a generation is better than another one if its evaluation is greater, ties being broken by index:

def better (j k : Fin lam) : Prop := evals k < evals j ∨ (evals j = evals k ∧ j < k)

Rather than sorting the generation, which would require to manipulate a permutation of \{1, \dots, \lambda\}, we count, for each point, the number of points that are better than it:

noncomputable def rank (k : Fin lam) : ℕ := #{j | better evals j k}

def

CMAES.rank {lam : ℕ} (evals : Fin lam → ℝ) (k : Fin lam) : ℕ
CMAES.rank {lam : ℕ} (evals : Fin lam → ℝ)
  (k : Fin lam) : ℕ

The rank of the k-th point of a generation, i.e. the number of points of that generation that are CMAES.better than it.

As better is a strict total order, this is a bijection between the points of the generation and \{0, \dots, \lambda - 1\}: the rank of a point is its index in the sorted generation, the best point having rank 0.

2.5.1.2. The state and the strategy parameters🔗

def

CMAES.State (d : ℕ) : Type
CMAES.State (d : ℕ) : Type

The state of CMA-ES: the mean m, the step size σ, the covariance matrix C and the evolution paths p_c and p_σ. The associated sampling distribution is 𝓝(m, σ² C).

structure

CMAES.Params : Type
CMAES.Params : Type

The strategy parameters of CMA-ES, i.e. the constants it does not adapt.

Constructor

CMAES.Params.mk

Fields

w : ℕ → ℝ

The recombination weights: w i is the weight of the point of rank i. They usually sum to one for i < μ and vanish for i ≥ μ, so that only the μ best points are recombined.

c_σ : ℝ

The learning rate of p_σ.

d_σ : ℝ

The damping of the step size update.

c_c : ℝ

The learning rate of p_c.

c_1 : ℝ

The learning rate of the rank-one update of C.

c_μ : ℝ

The learning rate of the rank-μ update of C.

The strategy parameters are constants: only the state is adapted along the iterations. The usual values, which depend on the dimension and on the size of the generations, are given by:

def

CMAES.defaultParams (d lam : ℕ) : CMAES.Params
CMAES.defaultParams (d lam : ℕ) :
  CMAES.Params

The default strategy parameters in dimension d for a generation of λ points, as suggested in (The CMA Evolution Strategy: A Tutorial, Hansen, 2023).

2.5.1.3. Updating the state🔗

Writing x_{i:\lambda} for the point of rank i, y_{i:\lambda} = (x_{i:\lambda} - m) / \sigma for its step and \langle y \rangle_w = \sum_{i = 1}^{\mu} w_i y_{i:\lambda} for the weighted recombination of the steps, the state is updated as \begin{aligned} m' &= m + \sigma \langle y \rangle_w, \\ p_\sigma' &= (1 - c_\sigma) p_\sigma + \sqrt{c_\sigma (2 - c_\sigma) \mu_{\text{eff}}} \; C^{-\frac{1}{2}} \langle y \rangle_w, \\ \sigma' &= \sigma \exp\left(\frac{c_\sigma}{d_\sigma} \left(\frac{\|p_\sigma'\|}{\mathbb{E}\|\mathcal{N}(0, I)\|} - 1\right)\right), \\ p_c' &= (1 - c_c) p_c + h_\sigma \sqrt{c_c (2 - c_c) \mu_{\text{eff}}} \; \langle y \rangle_w, \\ C' &= \left(1 - c_1 - c_\mu \sum_i w_i\right) C + c_1 \left(p_c' {p_c'}^\top + (1 - h_\sigma) c_c (2 - c_c) C\right) + c_\mu \sum_{i = 1}^{\mu} w_i y_{i:\lambda} y_{i:\lambda}^\top. \end{aligned} The new mean is the weighted recombination of the \mu best points, the step size is adapted along the conjugate evolution path p_\sigma (the cumulative step size adaptation), and the covariance matrix is the sum of a rank-one update, driven by the evolution path p_c, and of a rank-\mu update, driven by the steps of the generation.

Since the weights vanish beyond the \mu-th one, the sums over the sorted generation are simply sums over the generation, each point being weighted according to its rank:

noncomputable def weightedStep : EuclideanSpace ℝ (Fin d) := ∑ k, p.w (rank evals k) • step s pop k

One iteration gathers the five rules above:

noncomputable def update : State d :=
  (nextMean p s pop evals, nextStepSize p s pop evals, nextCov p g s pop evals,
    nextPathC p g s pop evals, nextPathσ p s pop evals)

def

CMAES.update {d lam : ℕ} (p : CMAES.Params) (g : ℕ) (s : CMAES.State d)
  (pop : ℝ_ d lam) (evals : Fin lam → ℝ) : CMAES.State d
CMAES.update {d lam : ℕ}
  (p : CMAES.Params) (g : ℕ)
  (s : CMAES.State d) (pop : ℝ_ d lam)
  (evals : Fin lam → ℝ) : CMAES.State d

One iteration of CMA-ES, from the population of the generation g and its evaluations.

2.5.1.4. The algorithm🔗

The state is a deterministic function of the past generations and of their evaluations, so that it can be recovered by recursion over the history of the algorithm:

noncomputable def state : (n : ℕ) → prod_iter_image (ℝ_ d lam) (Fin lam → ℝ) n → State d
  | 0, data => update p 0 s₀ (data.1 ⟨0, mem_Iic.mpr le_rfl⟩) (data.2 ⟨0, mem_Iic.mpr le_rfl⟩)
  | n + 1, data => update p (n + 1)
      (state n (Tuple.subTuple n.le_succ data.1, Tuple.subTuple n.le_succ data.2))
      (data.1 ⟨n + 1, mem_Iic.mpr le_rfl⟩) (data.2 ⟨n + 1, mem_Iic.mpr le_rfl⟩)

def

CMAES.state {d lam : ℕ} (p : CMAES.Params) (s₀ : CMAES.State d)
  (n : ℕ) : prod_iter_image (ℝ_ d lam) (Fin lam → ℝ) n → CMAES.State d
CMAES.state {d lam : ℕ} (p : CMAES.Params)
  (s₀ : CMAES.State d) (n : ℕ) :
  prod_iter_image (ℝ_ d lam) (Fin lam → ℝ)
      n →
    CMAES.State d

The state of CMA-ES once the generations 0, …, n have been sampled and evaluated, starting from s₀, i.e. the state from which the generation n + 1 is sampled.

The mean and the covariance matrix of the generation n + 1 being measurable functions of that state, the original CMA-ES is an instance of the above scheme, the evaluation space being \mathbb{R}^\lambda:

noncomputable def CMA_ES_original (p : Params) (s₀ : State d) :
    Algorithm (ℝ_ d lam) (Fin lam → ℝ) :=
  CMA_ES d lam (measurable_mean p s₀) (measurable_covar p s₀) s₀.m (s₀.σ ^ 2 • s₀.C)

def

CMA_ES_original {d lam : ℕ} (p : CMAES.Params) (s₀ : CMAES.State d) :
  Algorithm (ℝ_ d lam) (Fin lam → ℝ)
CMA_ES_original {d lam : ℕ}
  (p : CMAES.Params)
  (s₀ : CMAES.State d) :
  Algorithm (ℝ_ d lam) (Fin lam → ℝ)

The original CMA-ES algorithm for global optimization, starting from the state s₀ and using the strategy parameters p, e.g. CMAES.defaultParams. The λ points of each generation are sampled i.i.d. according to 𝓝(m, σ² C), the state being updated by CMAES.update. It is meant to be used with an evaluation function of the form fun x i ↦ f (x i), f being the objective function.