From 1735660ffe9d61098771048febb978d00c42417a Mon Sep 17 00:00:00 2001 From: Johann Dreo Date: Fri, 13 Jul 2012 14:58:27 +0200 Subject: [PATCH] use rank mu selector ; bugfix estimator's linear algebra : mu is useless in estimator ; arx = pop^T ; store D as a diagonal ; cwise prod for covar recomposition ; more asserts --- edo/application/cmaes/main.cpp | 6 ++-- edo/src/edoEstimatorNormalAdaptive.h | 52 +++++++++++++++++----------- edo/src/edoNormalAdaptive.h | 30 ++++++++-------- edo/src/edoSamplerNormalAdaptive.h | 20 ++++++----- 4 files changed, 63 insertions(+), 45 deletions(-) diff --git a/edo/application/cmaes/main.cpp b/edo/application/cmaes/main.cpp index 7e73d49dc..7a0e73f6d 100644 --- a/edo/application/cmaes/main.cpp +++ b/edo/application/cmaes/main.cpp @@ -57,7 +57,7 @@ int main(int ac, char** av) eoState state; // Instantiate all needed parameters for EDA algorithm - double selection_rate = parser.createParam((double)0.5, "selection_rate", "Selection Rate", 'R', section).value(); // R + //double selection_rate = parser.createParam((double)0.5, "selection_rate", "Selection Rate", 'R', section).value(); // R unsigned long max_eval = parser.getORcreateParam((unsigned long)0, "maxEval", "Maximum number of evaluations (0 = none)", 'E', "Stopping criterion").value(); // E @@ -69,10 +69,10 @@ int main(int ac, char** av) edoNormalAdaptive distribution(dim); - eoSelect< RealVec >* selector = new eoDetSelect< RealVec >( selection_rate ); + eoSelect< RealVec >* selector = new eoRankMuSelect< RealVec >( mu ); state.storeFunctor(selector); - edoEstimator< Distrib >* estimator = new edoEstimatorNormalAdaptive( distribution, mu ); + edoEstimator< Distrib >* estimator = new edoEstimatorNormalAdaptive( distribution ); state.storeFunctor(estimator); eoEvalFunc< RealVec >* plainEval = new Rosenbrock< RealVec >(); diff --git a/edo/src/edoEstimatorNormalAdaptive.h b/edo/src/edoEstimatorNormalAdaptive.h index 0e5a7b028..a2566772b 100644 --- a/edo/src/edoEstimatorNormalAdaptive.h +++ b/edo/src/edoEstimatorNormalAdaptive.h @@ -45,12 +45,11 @@ class edoEstimatorNormalAdaptive : public edoEstimatorAdaptive< EOD > { public: typedef typename EOT::AtomType AtomType; - typedef typename EOD::Vector Vector; + typedef typename EOD::Vector Vector; // column vectors @see edoNormalAdaptive typedef typename EOD::Matrix Matrix; - edoEstimatorNormalAdaptive( EOD& distrib, unsigned int mu ) : + edoEstimatorNormalAdaptive( EOD& distrib ) : edoEstimatorAdaptive( distrib ), - _mu(mu), _calls(0), _eigeneval(0) {} @@ -61,14 +60,15 @@ private: // compute recombination weights Eigen::VectorXd weights( pop_size ); double sum_w = 0; - for( unsigned int i = 0; i < _mu; ++i ) { - double w_i = log( _mu + 0.5 ) - log( i + 1 ); + for( unsigned int i = 0; i < pop_size; ++i ) { + double w_i = log( pop_size + 0.5 ) - log( i + 1 ); weights(i) = w_i; sum_w += w_i; } // normalization of weights weights /= sum_w; + assert( weights.size() == pop_size); return weights; } @@ -97,17 +97,18 @@ public: // Here, if we are in canonical CMA-ES, // pop is supposed to be the mu ranked better solutions, // as the rank mu selection is supposed to have occured. - Matrix arx( pop.size(), N ); + Matrix arx( N, lambda ); // copy the pop (most probably a vector of vectors) in a Eigen3 matrix - for( unsigned int i = 0; i < lambda; ++i ) { - for( unsigned int d = 0; d < N; ++d ) { - arx(i,d) = pop[i][d]; + for( unsigned int d = 0; d < N; ++d ) { + for( unsigned int i = 0; i < lambda; ++i ) { + arx(d,i) = pop[i][d]; // NOTE: pop = arx.transpose() } // dimensions } // individuals // muXone array for weighted recombination - Eigen::VectorXd weights = edoCMAESweights( N ); + Eigen::VectorXd weights = edoCMAESweights( lambda ); + assert( weights.size() == lambda ); // FIXME exposer les constantes dans l'interface @@ -137,6 +138,8 @@ public: Matrix invsqrtC = d.coord_sys() * d.scaling().asDiagonal().inverse() * d.coord_sys().transpose(); + assert( invsqrtC.innerSize() == d.coord_sys().innerSize() ); + assert( invsqrtC.outerSize() == d.coord_sys().outerSize() ); // expectation of ||N(0,I)|| == norm(randn(N,1)) double chiN = sqrt(N)*(1-1/(4*N)+1/(21*pow(N,2))); @@ -148,8 +151,11 @@ public: // compute weighted mean into xmean Vector xold = d.mean(); - d.mean( arx * weights ); - Vector xmean = d.mean(); + assert( xold.size() == N ); + // xmean ( N, 1 ) = arx( N, lambda ) * weights( lambda, 1 ) + Vector xmean = arx * weights; + assert( xmean.size() == N ); + d.mean( xmean ); /********************************************************************** @@ -157,9 +163,9 @@ public: *********************************************************************/ // cumulation for sigma - d.path_sigma( - (1.0-cs)*d.path_sigma() + sqrt(cs*(2.0-cs)*mueff)*invsqrtC*(xmean-xold)/d.sigma() - ); + d.path_sigma( + (1.0-cs)*d.path_sigma() + sqrt(cs*(2.0-cs)*mueff)*invsqrtC*(xmean-xold)/d.sigma() + ); // sign of h double hsig; @@ -173,10 +179,16 @@ public: // cumulation for the covariance matrix d.path_covar( - (1.0-cc)*d.path_covar() + hsig*sqrt(cc*(2.0-cc)*mueff)*(xmean-xold) / d.sigma() - ); + (1.0-cc)*d.path_covar() + hsig*sqrt(cc*(2.0-cc)*mueff)*(xmean-xold) / d.sigma() + ); - Matrix artmp = (1.0/d.sigma()) * arx - xold.rowwise().replicate(_mu); + Matrix xmu( N, lambda); + xmu = xold.rowwise().replicate(lambda); + assert( xmu.innerSize() == N ); + assert( xmu.outerSize() == lambda ); + Matrix artmp = (1.0/d.sigma()) * (arx - xmu); + // Matrix artmp = (1.0/d.sigma()) * arx - xold.colwise().replicate(lambda); + assert( artmp.innerSize() == N && artmp.outerSize() == lambda ); /********************************************************************** @@ -214,10 +226,11 @@ public: Eigen::SelfAdjointEigenSolver eigensolver( d.covar() ); // FIXME use JacobiSVD? d.coord_sys( eigensolver.eigenvectors() ); Matrix D = eigensolver.eigenvalues().asDiagonal(); + assert( D.innerSize() == N && D.outerSize() == N ); // from variance to standard deviations D.cwiseSqrt(); - d.scaling( D ); + d.scaling( D.diagonal() ); } return d; @@ -225,7 +238,6 @@ public: protected: - unsigned int _mu; unsigned int _calls; unsigned int _eigeneval; diff --git a/edo/src/edoNormalAdaptive.h b/edo/src/edoNormalAdaptive.h index 9ceb6b9fc..85db27b53 100644 --- a/edo/src/edoNormalAdaptive.h +++ b/edo/src/edoNormalAdaptive.h @@ -41,10 +41,11 @@ class edoNormalAdaptive : public edoDistrib< EOT > public: //typedef EOT EOType; typedef typename EOT::AtomType AtomType; - typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; // column vectors ( n lines, 1 column) typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix; edoNormalAdaptive( unsigned int dim ) : + _dim(dim), _mean( Vector::Zero(dim) ), _C( Matrix::Identity(dim,dim) ), _B( Matrix::Identity(dim,dim) ), @@ -53,7 +54,7 @@ public: _p_c( Vector::Zero(dim) ), _p_s( Vector::Zero(dim) ) { - assert( dim > 0); + assert( _dim > 0); } edoNormalAdaptive( unsigned int dim, @@ -88,23 +89,24 @@ public: return _mean.innerSize(); } - Vector mean() const {return _mean;} - Matrix covar() const {return _C;} - Matrix coord_sys() const {return _B;} - Vector scaling() const {return _D;} - double sigma() const {return _sigma;} + Vector mean() const {return _mean;} + Matrix covar() const {return _C;} + Matrix coord_sys() const {return _B;} + Vector scaling() const {return _D;} + double sigma() const {return _sigma;} Vector path_covar() const {return _p_c;} Vector path_sigma() const {return _p_s;} - void mean( Vector m ) { _mean = m; } - void covar( Matrix c ) { _C = c; } - void coord_sys( Matrix b ) { _B = b; } - void scaling( Vector d ) { _D = d; } - void sigma( double s ) { _sigma = s; } - void path_covar( Vector p ) { _p_c = p; } - void path_sigma( Vector p ) { _p_s = p; } + void mean( Vector m ) { _mean = m; assert( m.size() == _dim ); } + void covar( Matrix c ) { _C = c; assert( c.innerSize() == _dim && c.outerSize() == _dim ); } + void coord_sys( Matrix b ) { _B = b; assert( b.innerSize() == _dim && b.outerSize() == _dim ); } + void scaling( Vector d ) { _D = d; assert( d.size() == _dim ); } + void sigma( double s ) { _sigma = s; assert( s != 0.0 );} + void path_covar( Vector p ) { _p_c = p; assert( p.size() == _dim ); } + void path_sigma( Vector p ) { _p_s = p; assert( p.size() == _dim ); } private: + unsigned int _dim; Vector _mean; // Matrix _C; // covariance matrix Matrix _B; // eigen vectors / coordinates system diff --git a/edo/src/edoSamplerNormalAdaptive.h b/edo/src/edoSamplerNormalAdaptive.h index edf531eff..e8188105a 100644 --- a/edo/src/edoSamplerNormalAdaptive.h +++ b/edo/src/edoSamplerNormalAdaptive.h @@ -58,24 +58,28 @@ public: EOT sample( EOD& distrib ) { - unsigned int size = distrib.size(); - assert(size > 0); + unsigned int N = distrib.size(); + assert( N > 0); // T = vector of size elements drawn in N(0,1) - Vector T( size ); - for ( unsigned int i = 0; i < size; ++i ) { + Vector T( N ); + for ( unsigned int i = 0; i < N; ++i ) { T( i ) = rng.normal(); } - assert(T.innerSize() == size); + assert(T.innerSize() == N ); assert(T.outerSize() == 1); - Vector sol = distrib.mean() + distrib.sigma() * distrib.coord_sys() * (distrib.scaling().dot(T) ); + // mean(N,1) + sigma * B(N,N) * ( D(N,1) .* T(N,1) ) + Vector sol = distrib.mean() + + distrib.sigma() + * distrib.coord_sys() * (distrib.scaling().cwiseProduct(T) ); // C * T = B * (D .* T) + assert( sol.size() == N ); /*Vector sol = distrib.mean() + distrib.sigma() * distrib.coord_sys().dot( distrib.scaling().dot( T ) );*/ // copy in the EOT structure (more probably a vector) - EOT solution( size ); - for( unsigned int i = 0; i < size; i++ ) { + EOT solution( N ); + for( unsigned int i = 0; i < N; i++ ) { solution[i]= sol(i); }