CMASolutions reports best candidate's x value incorrect

[godrays]

Hello,

I found that the CMASolutions does not report best candidate's x value properly.

Here is an example output from 20 iteration:

Best solution: f-value: -6, x: 8.94434
Best solution: f-value: -6, x: 8.41353
Best solution: f-value: -6, x: 8.41353
Best solution: f-value: -6, x: 8.12711
Best solution: f-value: -6, x: 8.12711
Best solution: f-value: -6, x: 8.94491
Best solution: f-value: -6, x: 8.94491
Best solution: f-value: -6, x: 8.56085
Best solution: f-value: -6, x: 8.56085
Best solution: f-value: -6, x: 8.01153
Best solution: f-value: -6, x: 8.21036
Best solution: f-value: -6, x: 8.21036
Best solution: f-value: -6, x: 37.9178
Best solution: f-value: -6, x: 37.9178
Best solution: f-value: -6, x: 37.1124
Best solution: f-value: -6, x: 8.26821
Best solution: f-value: -6, x: 8.26821
Best solution: f-value: -6, x: 8.59635
Best solution: f-value: -6, x: 8.59635
Best solution: f-value: -6, x: 8.84206

I validated that x values in fitness function never gets higher than 21 during simulations, which is expected. You can see some of the x values are over 21. Sometimes I see negative x values too, which is also NOT expected.

Here is the test code I used:

#include <iostream>
#include <limits>
#include <libcmaes/cmaes.h>

using namespace libcmaes;

int main()
{
                      // x: 0  1  2  3   4   5   6   7   8   9  10 11 12 13 14 15 16 17  18  19  20, 21
    std::vector<double>  y{ 3, 2, 1, 0, -1, -2, -3, -4, -6, -3, -1, 2, 5, 7, 8, 9, 8, 3, -2, -4, -5, -1 };

    int minX = std::numeric_limits<int>::max();
    int maxX = std::numeric_limits<int>::min();

    FitFunc fsphere = [&](const double *x, const int N)
    {
        int i = static_cast<int>(x[0]);

        minX = std::min<int>(i, minX);
        maxX = std::max<int>(i, maxX);

        return y[i];
    };

    const int dim = 1;
    const double sigma = 1;
    const int numOfRun = 20;    // run the same optimization n times

    double lbounds[dim];        // lower parameter bounds
    double ubounds[dim];        // upper parameter bounds

    for (int n=0; n<numOfRun; ++n)
    {
        std::vector<double>  x0(dim, 15.0);

        for (int i=0; i<dim; ++i)
        {
            lbounds[i] = 0.0;
            ubounds[i] = static_cast<double>(y.size());
        }

        GenoPheno<pwqBoundStrategy> gp(lbounds, ubounds, dim);
        CMAParameters<GenoPheno<pwqBoundStrategy>> cmaparams(x0, sigma, 100, 0, gp);
        cmaparams.set_algo(aCMAES);
        CMASolutions cmasols = cmaes<GenoPheno<pwqBoundStrategy>>(fsphere, cmaparams);

        std::cout << "Best solution: f-value: " << cmasols.best_candidate().get_fvalue()
                  << ", x: " << cmasols.best_candidate().get_x().front() << std::endl;
    }

    std::cout << "min X: " << minX << ", max X: " << maxX << std::endl;

    return 0;
}

Is this expected behavior? If yes then why?

I compiled the lib with OpenMP=OFF
I use v0.10 and built it from the source code.

Thanks

[godrays]

I found a suggestion in other issues saying we needed to use cmasols.get_best_seen_candidate() instead of cmasols.best_candidate()

    auto best_candidate = cmasols.get_best_seen_candidate();
    
    std::cout << "Best solution: f-value: " << best_candidate.get_fvalue()
              << ", x: " << best_candidate.get_x().front() << std::endl;

Unfortunately, I still see incorrect answers in outputs.. Just see below that x should never be lower than 0.

...
Best solution: f-value: -6, x: 8.78817
Best solution: f-value: -6, x: -8.36663

Here is an other incorrect output.. x should never be higher than 21.

...
Best solution: f-value: -6, x: 8.6857
Best solution: f-value: -6, x: 8.69536
Best solution: f-value: -6, x: 37.5472
Best solution: f-value: -6, x: 37.5472
Best solution: f-value: -6, x: 8.07772
...

Please note that the issue happened less frequently than the first test, but still the same issue.

[beniz]

Hi, see https://github.com/CMA-ES/libcmaes/wiki/Defining-and-using-bounds-on-parameters
You are using a genopheno for bounded parameters, thus need to print the solution in pheno space: cmasols.print(std::cout, 0, gp) ?
Also, for you own better understanding, the implementation is here: https://github.com/CMA-ES/libcmaes/blob/master/src/cmasolutions.cc#L257 (see the gp.get_pheno(best_candidate()... call.

[godrays]

Thank you for the quick clarification.

Transformation solved my issue. However, cmasols.print(...) implementation you mentioned uses the best_candidate() for report

<< " / x=" << gp.pheno(best_candidate().get_x_dvec()).transpose();

The documentation suggests using best_seen_candidate(), which is confusing.

Eigen::VectorXd bestparameters = gp.pheno(cmasols.get_best_seen_candidate().get_x_dvec());

In this case, anyone who uses cmasols.print(...) might not see the best parameters. What is the exact difference between best_candidate() and best_seen_candidate()? Can we document the diff?

[beniz]

best_seen_candidate() returns the best over the whole run, while best_candidate() returns the current solution.

You should be able to apply gp.pheno to best_seen_candidate(), this is the reason why I pointed to the implementation.

Maybe we should provide a print_best_seen function as well...

[mirai-computing]

Hi! Same problem here.
Both best_candidate() and best_seen_candidate() return solution vector via get_x() in some weird local scale, so I had to resort to selecting best solution at every evaluation of the objective function parallel to libcmaes inner workings.
The problem setup is following:

libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy,libcmaes::linScalingStrategy>> cmaparams(x0,sigmas);
cmaparams.set_algo(aCMAES);
libcmaes::CMASolutions cmasols = libcmaes::cmaes<libcmaes::GenoPheno<libcmaes::NoBoundStrategy,libcmaes::linScalingStrategy>>(f,cmaparams);
return cmasols.best_seen_candidate().get_x();

How do I get correct gp.pheno() transformation applied to result in this setup?
In print function gp is private argument
template <class TGenoPheno> std::ostream& CMASolutions::print(std::ostream &out, const int &verb_level, const TGenoPheno &gp) const
so doing this transformation on outside doesn't seem possible
" / x=" << gp.pheno(best_candidate().get_x_dvec()).transpose()

It would be nice to provide std::vector<double> get_x_global() to Candidate that would return correctly transformed version of solution vector so that it would contain exactly same values that were used to evaluate objective function.

[beniz]

Hello @mirai-computing the gp (GenoPheno<NoBoundStrategy>) should be defined in user-space, have you looked at this wiki page ? https://github.com/CMA-ES/libcmaes/wiki/Using-parameter-space-transforms-known-as-genotype-phenotype-transforms

[mirai-computing]

Thanks, this clears things up a bit.
I haven't seen this page yet since I'm new to libcmaes and just trying to understand how everything works and heavy use of C++ templates makes it hard to figure out what should go where by design etc. I'm kind of old style person preferring runtime switches and pointers over compile time specifications :)

[mirai-computing]

Now there's another problem. When I instantiate gp and parameters like

libcmaes::GenoPheno<libcmaes::NoBoundStrategy,libcmaes::linScalingStrategy> gp;
libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy,libcmaes::linScalingStrategy>> cmaparams(x0,sigma,-1,0,gp);

with sigma being a single number and not a vector of same size as problem dimension,
it starts failing at me with with following message as soon as dimension is 2.

/usr/include/eigen3/Eigen/src/Core/CwiseBinaryOp.h:116: Eigen::CwiseBinaryOp<BinaryOp, Lhs, Rhs>::CwiseBinaryOp(const Lhs&, const Rhs&, const BinaryOp&) [with BinaryOp = Eigen::internal::scalar_difference_op<double, double>; LhsType = const Eigen::Matrix<double, -1, 1>; RhsType = const Eigen::Matrix<double, -1, 1>; Eigen::CwiseBinaryOp<BinaryOp, Lhs, Rhs>::Lhs = Eigen::Matrix<double, -1, 1>; Eigen::CwiseBinaryOp<BinaryOp, Lhs, Rhs>::Rhs = Eigen::Matrix<double, -1, 1>]: Assertion `aLhs.rows() == aRhs.rows() && aLhs.cols() == aRhs.cols()' failed.

I've came across this issue earlier and I thought I've been calling wrong constructor since switching to one that allowed for anisotropic distribution with different sigmas (that were equal but came in a vector) fixed the issue.
But there's no CMAParameters constructor to pass both gp and std::vector<double>& sigmas at the same time.
What am I doing wrong?

Looking at CMAParamters sources I see there's more ways to pass parameters to the upper class but I'm not sure if this configuration I need (with both vector sigma and gp) makes sense to CMA-ES.

[mirai-computing]

Update: here's the call stack inside libcmaes and Eigen when the assert fails:

#5  0x00000000006142c7 in Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>::CwiseBinaryOp (func=..., aRhs=..., aLhs=..., this=<optimized out>) at /usr/include/eigen3/Eigen/src/Core/CwiseBinaryOp.h:116
#6  Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::operator-<Eigen::Matrix<double, -1, 1, 0, -1, 1> > (other=..., this=<optimized out>) at /usr/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseBinaryOps.h:19
#7  libcmaes::linScalingStrategy::scale_to_f (y=..., x=..., this=0x7fffffff9a90) at ../include/libcmaes/scaling.h:131
#8  libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy>::pheno (this=this@entry=0x7fffffff9a40, candidates=...) at ../include/libcmaes/genopheno.h:330
#9  0x0000000000617c58 in libcmaes::CMAStrategy<libcmaes::ACovarianceUpdate, libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> >::optimize(std::function<void (Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&)> const&, std::function<Eigen::Matrix<double, -1, -1, 0, -1, -1> ()> const&, std::function<void ()> const&) (this=0x7fffffff96c0, evalf=..., askf=..., tellf=...) at cmastrategy.cc:356
#10 0x00000000005c78c4 in libcmaes::CMAStrategy<libcmaes::ACovarianceUpdate, libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> >::optimize (this=0x7fffffff96c0) at lib/libcmaes/cmastrategy.h:114
#11 0x00000000005c5e2d in libcmaes::ESOptimizer<libcmaes::CMAStrategy<libcmaes::ACovarianceUpdate, libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> >, libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> >, libcmaes::CMASolutions>::optimize (this=0x7fffffff96c0) at lib/libcmaes/esoptimizer.h:115
#12 0x00000000005c3301 in libcmaes::cmaes<libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> >(std::function<double (double const*, int const&)>&, libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> >&, std::function<int (libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> > const&, libcmaes::CMASolutions const&)>&, std::function<Eigen::Matrix<double, -1, 1, 0, -1, 1> (double const*, int const&)>, libcmaes::CMASolutions const&, std::function<int (libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy, libcmaes::linScalingStrategy> > const&, libcmaes::CMASolutions const&, std::basic_ofstream<char, std::char_traits<char> >&)>&) (func=..., parameters=..., pfunc=..., gfunc=..., solutions=..., pffunc=...) at lib/libcmaes/cmaes.h:123

[mirai-computing]

This looks like a problem at the interface with Eigen library.
The problem begins in libcmaes/genopheno.h when it calls scale_to_f():

      for (int i=0;i<ncandidates.cols();i++)
	{
	  dVec ycoli;
	  _scalingstrategy.scale_to_f(ncandidates.col(i),ycoli);
	  ncandidates.col(i) = ycoli;
	}

Then linScalingStrategy (libcmaes/scaling.h) tries to subtract two vectors x and _shift:

    void scale_to_f(const dVec &x,
		    dVec &y) const
    {
      y = x - _shift;
      y = y.cwiseQuotient(_scaling);
    }

And here Eigen fails since for some reason x and _shift vectors have different sizes

    CwiseBinaryOp(const Lhs& aLhs, const Rhs& aRhs, const BinaryOp& func = BinaryOp())
      : m_lhs(aLhs), m_rhs(aRhs), m_functor(func)
    {
      EIGEN_CHECK_BINARY_COMPATIBILIY(BinaryOp,typename Lhs::Scalar,typename Rhs::Scalar);
      // require the sizes to match
      EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
      eigen_assert(aLhs.rows() == aRhs.rows() && aLhs.cols() == aRhs.cols());
    }

Any ideas?

[mirai-computing]

Update: debug print shows that _shift is wrong size: scale_to_f: y = x[1,2] - _shift[1,1].

[beniz]

libcmaes::CMAParameters<libcmaes::GenoPhenolibcmaes::NoBoundStrategy,libcmaes::linScalingStrategy> cmaparams(x0,sigma,-1,0,gp);

Aren't you missing the dim parameter to cmaparams ?

https://github.com/CMA-ES/libcmaes/wiki/Defining-and-using-bounds-on-parameters
https://github.com/CMA-ES/libcmaes/blob/master/src/cmaparameters.cc#L29

[mirai-computing]

Aren't you missing the dim parameter to cmaparams ?

Since I'm passing x0 as std::vector<double> with a certain size and not a pointer to array of doubles, I expect cmaparams to extract dim from x0.size() by itself. Is it wrong to do it this way? Should I only pass a pair of parameters int dim, double* x0.front() if I want it to explicitly know the dimension?

Then I wonder how is the dimension of problem determined when dim is not even in list of arguments? Like, 2 of 3 constructors of CMAParameters don't even have dim argument but accept vectors for x0 and sigma...
http://beniz.github.io/libcmaes/doc/html/classlibcmaes_1_1CMAParameters.html

[mirai-computing]

And no, this doesn't help. I switched call to libcmaes::CMAParameters<libcmaes::GenoPheno<libcmaes::NoBoundStrategy,libcmaes::linScalingStrategy>> cmaparams(dim,x0ptr,sigma,-1,0,gp); syntax and it still fails when dim>1 since size of _shift is still 1.

[mirai-computing]

p.s. this worked perfectly without NoBoundStrategy,linScalingStrategy specification, but if I get back to default spec then I loose scale adaptation capabilities of CMA-ES which is bad in its own right. But turning this spec on somehow makes libcmaes to loose dimensional awareness. I've only found that without bounds _shift initializes to 1-vector in constructor of linScalingStrategy and it seems it is not reallocated to dim size later.

[mirai-computing]

OK, I kind of solved it.
Turns out it was necessary for GenoPheno with NoBoundStrategy,linScalingStrategy spec to get to a constructor with externally supplied scaling and shift dim-size vectors and declared but unused lbounds,ubounds arguments.

 dVec scaling = dVec::Constant(dim,1.0);
 dVec shift = dVec::Zero(dim);
 // lbounds,unbounds are unused with NoBoundStrategy, so we can call gp() with fake null pointers.
 libcmaes::GenoPheno<libcmaes::NoBoundStrategy,libcmaes::linScalingStrategy> gp(scaling,shift,0,0);

Since there's no default GenoPheno constructor for NoBoundStrategy,linScalingStrategy spec, it ends up in a default constructor for NoBoundStrategy,NoScalingStrategy spec that only initializes _scaling and _shift to 1-vectors. And this leads to mismatch between selected scaling strategy and internal state which results is supplying incorrectly sized vectors to Eigen where it triggers assert violation and fails.

[beniz]

Great, perseverance always pays :)
I haven´t looked at it in a long time, but the constructor should be here: https://github.com/CMA-ES/libcmaes/blob/master/include/libcmaes/genopheno.h#L336
If something is wrong, feel free to PR a fix / addition.

"I'm kind of old style person preferring runtime switches and pointers over compile time specifications :)"

I'm not too young either :) When I wrote libcmaes I was pushing policy-based design as much as I could. I don't regret it, and this type of design, that avoids many pointer stuff and errors at runtime, has stayed with me until today. However, I did push it a bit too much in libcmaes, making the C++ API a bit difficult to use.

That being said, it's been pretty stable over years, and I like having a function optimizer running on safe grounds, controled memory, etc.. as offered by policy-based templates!

[mirai-computing]

Well, you're cool and libcmaes is great :)

I do understand why it is created in this way and I appreciate it, my only complaint up till now is that it is hard to debug things when almost everything important happens during static initialization in constructors ;>.<;
However, seriously speaking, there's at least one case where this style may encounter a huge problem and it is purely technical issue and it is time and execution flow control. For a task that may consume a significant amount of CPU time unpredictably, it is beneficial to have a state machine that does one short time limited action at a time. Multithreading an app partially does this from CPU time consumption control side of things but it is bad at having a definitive state of a halted thread (not algorithmically synchronized to underlying computational method) :)

I'm trying to marry CMA-ES in shape of libcmaes implementation to a SR (symbolic regression) implementation as a more of less universal solver for free parameters and I must say it converges pretty well and fast at least in not too complicated cases. So, I confirm you did a great job.