[WIP] add initial version of CMA-ES with Margin

[nomuramasahir0]

The initial implementation of CMA-ES with Margin. I haven't had enough time to test it carefully, but I have checked that it works properly for the script below.

import numpy as np
from cmaes import CMAwM


def ellipsoid_onemax(x, n_zdim):
    n = len(x)
    n_rdim = n - n_zdim
    r = 10
    if len(x) < 2:
        raise ValueError("dimension must be greater one")
    ellipsoid = sum([(1000 ** (i / (n_rdim - 1)) * x[i]) ** 2 for i in range(n_rdim)])
    onemax = n_zdim - (0.0 < x[(n - n_zdim) :]).sum()
    return ellipsoid + r * onemax


def main():
    dim = 20
    binary_dim = dim // 2
    discrete_space = np.tile(
        np.arange(0, 2, 1), (binary_dim, 1)
    )  # binary variables (dim_bi x 2)
    optimizer = CMAwM(mean=0 * np.ones(dim), sigma=2.0, discrete_space=discrete_space)
    print(" evals    f(x)")
    print("======  ==========")

    evals = 0
    while True:
        solutions = []
        for _ in range(optimizer.population_size):
            x_for_eval, x_for_tell = optimizer.ask()
            value = ellipsoid_onemax(x_for_eval, binary_dim)
            evals += 1
            solutions.append((x_for_tell, value))
            if evals % 300 == 0:
                print(f"{evals:5d}  {value:10.5f}")
        optimizer.tell(solutions)

        if optimizer.should_stop():
            break


if __name__ == "__main__":
    main()

[nomuramasahir0]

@c-bata I just published the first version of the CMA-ES with Margin implementation. I’d be happy if i could get your feedback (, though I may not have much time to address it).

[dietmarwo]

Very interesting work. Maybe I will integrate something similar in my own CMA implementation. Differential Evolution suffers from a similar issue, which I already fixed (see below). BiteOpt seems to work "out of the box". CR-FM-NES also has the problem. For comparison you may try:

from scipy.optimize import Bounds
from fcmaes.optimizer import De_cpp, Bite_cpp, Crfmnes_cpp, wrapper
from fcmaes import retry

# do 'pip install fcmaes'
    
def ellipsoid_onemax(x, n_zdim):
    n = len(x)
    n_rdim = n - n_zdim
    r = 10
    if len(x) < 2:
        raise ValueError("dimension must be greater one")
    ellipsoid = sum([(1000 ** (i / (n_rdim - 1)) * x[i]) ** 2 for i in range(n_rdim)])
    onemax = n_zdim - (0.0 < (x[(n - n_zdim) :]).astype(int)).sum()
    return ellipsoid + r * onemax

def main():
    
    dim = 20
    n_zdim = dim // 2

    lb = [-1000]*(dim-n_zdim) + [0]*n_zdim
    ub = [1000]*(dim-n_zdim) + [1.9999]*n_zdim
    bounds = Bounds(lb, ub)

    def fun(x):
        return ellipsoid_onemax(x, n_zdim)
    
    fit = wrapper(fun)
    #opt = Bite_cpp(10000) # https://github.com/avaneev/biteopt works quite well
    #opt = De_cpp(10000) # differential evolution has a similar issue
    opt = De_cpp(6000, ints=[False]*(dim-n_zdim) + [True]*n_zdim) # fixed using the "ints" parameter
    #opt = Crfmnes_cpp(6000) # also has the issue (not fixed yet)
    retry.minimize(fit, bounds, num_retries=1, optimizer=opt)

if __name__ == "__main__":
    main()

[nomuramasahir0]

@dietmarwo Thanks for interesting information! It is interesting to compare the CMA-ES with Margin with other handling methods including yours. I believe the idea of CMA-ES with Margin can be extended to CR-FM-NES, and it is work in progress.

[knshnb]

@nomuramasahir0 Hi, I appreciate your implementation and I'm looking forward to integrating this feature into Optuna! Does this PR currently support integer input? In that case, is the below usage correct, where CMAwM does not seem to work correctly as CMA?

from cmaes import CMA
from cmaes import CMAwM
import numpy as np


def sphere(x) -> float:
    return (x[0] - 2) ** 2 + round(x[1] - 2) ** 2


def main(optimizer):
    print(" evals    f(x)")
    print("======  ==========")

    evals = 0
    while True:
        solutions = []
        for _ in range(optimizer.population_size):
            if isinstance(optimizer, CMAwM):
                x_for_eval, x_for_tell = optimizer.ask()
            else:
                x_for_eval = optimizer.ask()
                x_for_tell = x_for_eval
            value = sphere(x_for_eval)
            evals += 1
            solutions.append((x_for_tell, value))
            if evals % 100 == 0:
                print(f"{evals:5d}  {value:10.5f}")
        optimizer.tell(solutions)

        if optimizer.should_stop():
            break


if __name__ == "__main__":
    optimizer = CMAwM(
        mean=np.zeros(2),
        sigma=2.0,
        discrete_space=np.array([[0, 10]]),
        continuous_space=np.array([[0.0, 10.0]]),
    )
    main(optimizer)
    optimizer = CMA(mean=np.zeros(2), sigma=2.0, bounds=np.array([[0.0, 10.0], [0.0, 10.0]]))
    main(optimizer)

 evals    f(x)
======  ==========
  100     4.00391
  200     4.00000
  300     4.00000
  400     4.00000
  500     4.00000
  600     4.00000
  700     4.00000
  800     4.00000
  900     4.00000
 1000     4.00000
 1100     4.00000
 1200     4.00000
 evals    f(x)
======  ==========
  100     0.00378
  200     0.00000
  300     0.00000
  400     0.00000

[nomuramasahir0]

@knshnb Thank you for having interest! discrete_space=np.array([[0, 10]]) means that the discrete space is {0, 10}. If you want to set the discrite space to {0, 1, 2, ..., 9, 10}, as well as the space in the CMA code, please rewrite the discrete_space in the CMAwM to

discrete_space=np.array([np.arange(0, 11, 1)])

But I must admit I haven't tested the implementation carefully, so there may be some bugs. I will investigate the validity of the implementation in the next two or three weeks.

[knshnb]

Thanks for the quick reply. Now it works correctly with your suggestion!

Let me ask another question. The current implementation does not seem to support discrete search spaces with different numbers of elements (please see the script below). Is it technically possible to support such cases?

from cmaes import CMAwM
import numpy as np


def sphere(x) -> float:
    return round(x[0] - 2) ** 2 + round(x[1] - 2) ** 2


def main(optimizer):
    print(" evals    f(x)")
    print("======  ==========")

    evals = 0
    while True:
        solutions = []
        for _ in range(optimizer.population_size):
            x_for_eval, x_for_tell = optimizer.ask()
            value = sphere(x_for_eval)
            evals += 1
            solutions.append((x_for_tell, value))
            if evals % 100 == 0:
                print(f"{evals:5d}  {value:10.5f}")
        optimizer.tell(solutions)

        if optimizer.should_stop():
            break


if __name__ == "__main__":
    optimizer = CMAwM(
        mean=np.zeros(2),
        sigma=2.0,
        discrete_space=np.array([np.arange(0, 11), np.arange(0, 6)], dtype=object),
    )
    main(optimizer)

File "venv/lib/python3.9/site-packages/cmaes/_cmawm.py", line 222, in __init__
    self.z_lim = (self.z_space[:, 1:] + self.z_space[:, :-1]) / 2
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

[nomuramasahir0]

@knshnb

The current implementation does not seem to support discrete search spaces with different numbers of elements (please see the script below). Is it technically possible to support such cases?

Thank you for the nice question. We support the case and it's possible. I believe the following code works for the case.

discrete_space = np.full((2, 11), np.nan)
discrete_space[0, :11] = np.arange(0, 11)
discrete_space[1, :6] = np.arange(0, 6)

By the way, round in the sphere function can be removed because x_for_eval is already encoded into discrete space.

[knshnb]

Thank you for the quick response again! I could successfully handle my case with your advice.

FYI: I implemented a prototypical integration into Optuna (link) and took a brief benchmark on HPO-Bench problems. CMA-ES with margin ("_CmaEsSampler_NopPruner_1") tended to have better performance than the original CMA-ES with many trials (200~800+).

[c-bata]

Let me close this PR since we had an internal discussion and decided to introduce CMA-ES with Margin with an alternative user interface at #121. Thank you again for your huge effort.