CMA-ES

[VolodymyrOrlov]

Implements CMA-ES algorithm, as described in this paper. Along with CMA-ES, this PR extends argmin-math with additional functions:

• argsort: https://numpy.org/doc/stable/reference/generated/numpy.argsort.html
• eigSym: https://www.rdocumentation.org/packages/MGMM/versions/0.3.1/topics/eigSym
• np.random.standard_normal: https://numpy.org/doc/stable/reference/random/generated/numpy.random.standard_normal.html
• np.outer: https://numpy.org/doc/stable/reference/generated/numpy.outer.html
• np.take: https://numpy.org/doc/stable/reference/generated/numpy.take.html

[stefan-k]

This is huge! Thanks a lot, not just for CMA-ES, but also for the massive extension of argmin-math, which will certainly come in handy in the future.
Since this is quite a large PR, I'll have to take some time to review it. I hope to be able to review it in the course of the next week.

[VolodymyrOrlov]

Sounds good, @stefan-k ! Thanks for taking a look at the PR!

[Armavica]

Thank you for all your work! I have been wanting to implement this for a long time but I am glad that someone managed to do it. I am very much looking forward to being able to use CMAES from within argmin!

I tried to increase the number of iterations from 100 to 1000 in the cmaes example and the search fails all the time because at some point all the fitnesses become NaN, do you know where this might come from?

I also tried computing the minimum of the 10D rosenbrock function with the same parameters, and the Shur transform is failing but I am not sure why.

[VolodymyrOrlov]

Thank you for giving it a try, @Armavica! What backend are you using, vec/ndarray or nalgebra? It is likely Eigendecomposition routine. For vec I took all the decomposition code from https://commons.apache.org/proper/commons-math/, but maybe their implementation is not numerically stable

[Armavica]

I just changed the parameters in cmaes.rs and ran it with cargo run --release --example cmaes, so I think it is vec unless I am wrong?

[VolodymyrOrlov]

Yes, it is vec. I'll see what I can do to make it more stable

[Armavica]

It looks like there is a similar problem with the nalgebra backend, observed if you change the maximum number of iterations from 100 to 10000 in test_nalgebra_solver?

[stefan-k]

I only have a few minor comments. Unfortunately I'm not a user of nalgebra therefore I'll have to trust you on that ;)

After addressing these comments and the problems @Armavica pointed out I think this is ready to go!

[VolodymyrOrlov]

I only have a few minor comments. Unfortunately I'm not a user of nalgebra therefore I'll have to trust you on that ;)

After addressing these comments and the problems @Armavica pointed out I think this is ready to go!

Thank you for reviewing my PR! I've found a couple of numerical issues that lead to problems pointed out by @Armavica! I hope to have some changes for you guys some time by the mid next week.

[VolodymyrOrlov]

I've fixed issues discovered by @Armavica . The vector-based implementation is still inferior to ndarray and nalgebra, but at least the method converges for multidimensional Rosenbrock and doesn't panic for when number of iterations exceeds 100

[stefan-k]

I've fixed issues discovered by @Armavica . The vector-based implementation is still inferior to ndarray and nalgebra, but at least the method converges for multidimensional Rosenbrock and doesn't panic for when number of iterations exceeds 100

Thanks a lot once again! :) Just out of curiosity, in which way is the vec-based implementation inferior to the others?

I only have a very few minor comments. One maybe larger one is the use of bulk_cost instead of cost which enables parallel computation of the cost function for the entire population. Other than that I think this is good to go. Before merging, could you squash your commits a bit please? It doesn't need to be a single commit but grouping them in a couple of commits would be nice.

[VolodymyrOrlov]

I've fixed issues discovered by @Armavica . The vector-based implementation is still inferior to ndarray and nalgebra, but at least the method converges for multidimensional Rosenbrock and doesn't panic for when number of iterations exceeds 100

Thanks a lot once again! :) Just out of curiosity, in which way is the vec-based implementation inferior to the others?

I only have a very few minor comments. One maybe larger one is the use of bulk_cost instead of cost which enables parallel computation of the cost function for the entire population. Other than that I think this is good to go. Before merging, could you squash your commits a bit please? It doesn't need to be a single commit but grouping them in a couple of commits would be nice.

The vector and nalgebra backend requires more iterations to converge because Eigenvector decomposition implementation is not as good a the one that comes with LAPACK

You will notice that number of modified files went up with my latest changes. This is due to commits squashing.

[stefan-k]

The vector and nalgebra backend requires more iterations to converge because Eigenvector decomposition implementation is not as good a the one that comes with LAPACK

Thanks for the explanation!

You will notice that number of modified files went up with my latest changes. This is due to commits squashing.

I think something went wrong during squashing. It seems that merging main into your branch (rather than rebasing) caused an unfortunate starting point for squashing. Most commits were squashed into one of mine, which led to them being assigned to my user. I tried to fix this, so now I'm kindly asking you to review the resulting changes once the CI is successful ;)
I saved the previous commits just in case I messed up.

[VolodymyrOrlov]

The vector and nalgebra backend requires more iterations to converge because Eigenvector decomposition implementation is not as good a the one that comes with LAPACK

Thanks for the explanation!

You will notice that number of modified files went up with my latest changes. This is due to commits squashing.

I think something went wrong during squashing. It seems that merging main into your branch (rather than rebasing) caused an unfortunate starting point for squashing. Most commits were squashed into one of mine, which led to them being assigned to my user. I tried to fix this, so now I'm kindly asking you to review the resulting changes once the CI is successful ;) I saved the previous commits just in case I messed up.

Thank you for your help with squashing commits, @stefan-k ! The changes looks good, feel free to merge this PR when you deem it ready.

[Armavica]

Hi, I pulled the changes from the branch VolodymirOrlov/main and modified the examples/cmaes.rs file to allow for 1000 iterations. The optimization still panics. Am I doing something wrong?

[stefan-k]

Hi, I pulled the changes from the branch VolodymirOrlov/main and modified the examples/cmaes.rs file to allow for 1000 iterations. The optimization still panics. Am I doing something wrong?

Interestingly, I have no issues with 1000 iterations, but I do see a panic for 2000 iterations:

thread 'main' panicked at 'Exceeded maximum number of iterations when calculating QL transformation', /home/stefan/coding/argmin/argmin/argmin-math/src/vec/eigen.rs:422:1

Edit: I'm unsure if it is sensible to increase EIGEN_MAX_ITERATIONS.

Edit 2: It is also not consistently reproducible, which makes sense given the used randomness.

[Armavica]

Ah yes indeed sorry, I was apparently doing something wrong: it also works most of the time for me now. I also get this panic sometimes, non-reproducibly. From a few rapid tests this implementation seems to work even even for high-dimensional functions. Congrats for fixing the issue!

The author of CMA-ES has a page with a few different examples to help test and benchmark CMA-ES implementations: https://cma-es.github.io/cmaes_sourcecode_page.html#testing I think it could be nice ideally to see how this implementation performs on these tests, but maybe this can come later…

[stefan-k]

From a quick look I've noticed that there are some direct comparisons of floats, for instance:

if (e[i + 1] == 0.0) {
    eigenvalues[i + 1] -= u;
    e[m] = 0.0;
    break;
}

and

if t == 0. && i >= j {
    continue;
}

Could this cause issues because it is unlikely that those will be exactly 0?

[stefan-k]

Hi @VolodymyrOrlov , this PR is really far ahead and I'd really like to include it in the next release. There only seem to be one remaining issue regarding the Eigenvalue decomposition for the vector and nalgebra backends and the PR would have to be rebased onto the current main. I'm unsure how to proceed. I wouldn't mind dropping support for the Vec backend as it is rather limited and inefficient anyways, but it would be sad if all the tremendous effort you put into this would be for nothing.
Support for the nalgebra backend however would be pretty great. What do you suggest?

[VolodymyrOrlov]

Hi @VolodymyrOrlov , this PR is really far ahead and I'd really like to include it in the next release. There only seem to be one remaining issue regarding the Eigenvalue decomposition for the vector and nalgebra backends and the PR would have to be rebased onto the current main. I'm unsure how to proceed. I wouldn't mind dropping support for the Vec backend as it is rather limited and inefficient anyways, but it would be sad if all the tremendous effort you put into this would be for nothing. Support for the nalgebra backend however would be pretty great. What do you suggest?

Hi @stefan-k! Thanks for continued support of this PR and my apologies for not acting sooner on suggested changes. I switched companies and now work for a small startup. This consumes most of my time and I didn't had a chance to do any open-source work lately. Give me some time tomorrow and early next week - I'd like to try to resolve remaining issues and update PR to prepare it for the merge.

[stefan-k]

Rereading my previous message again I realized that it may have come across as pushy. I apologize if it made that impression, this was certainly not my intention.
Thanks for continuing to work on this! Take your time, I know how difficult it can be to manage both work and open source contributions. Congratulations on the new job! :)

[Boscop]

@VolodymyrOrlov Thanks for your great work, I'm looking forward to this being merged :)