Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement Active-CMA-ES #235

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Implement Active-CMA-ES #235

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f682183
implement active updates in cma-es
egomeh Jun 17, 2018
df7c353
Make sure that active updates are default to false
egomeh Jun 17, 2018
ee22132
Add tests for active updates
egomeh Jun 17, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 30 additions & 0 deletions Test/Algorithms/DirectSearch/CMA.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,16 @@ BOOST_AUTO_TEST_CASE( CMA_Cigar )
testFunction( optimizer, function, 10, 1000, 1E-10 );
}

BOOST_AUTO_TEST_CASE( Active_CMA_Cigar )
{
Cigar function(4);
CMA optimizer;
optimizer.useActiveUpsates(true);

std::cout << "Testing: " << optimizer.name() << " with " << function.name() << " and active updates" << std::endl;
testFunction(optimizer, function, 10, 1000, 1E-10);
}

BOOST_AUTO_TEST_CASE( CMA_Discus )
{
Discus function(3);
Expand All @@ -81,6 +91,16 @@ BOOST_AUTO_TEST_CASE( CMA_Discus )
testFunction( optimizer, function, 10, 1000, 1E-10 );
}

BOOST_AUTO_TEST_CASE( Active_CMA_Discus )
{
Discus function(4);
CMA optimizer;
optimizer.useActiveUpsates(true);

std::cout << "\nTesting: " << optimizer.name() << " with " << function.name() << " and active updates" << std::endl;
testFunction(optimizer, function, 10, 1000, 1E-10);
}

BOOST_AUTO_TEST_CASE( CMA_Ellipsoid )
{
Ellipsoid function(5);
Expand All @@ -90,6 +110,16 @@ BOOST_AUTO_TEST_CASE( CMA_Ellipsoid )
testFunction( optimizer, function, 10, 1000, 1E-10 );
}

BOOST_AUTO_TEST_CASE( Active_CMA_Ellipsoid )
{
Ellipsoid function(10);
CMA optimizer;
optimizer.useActiveUpsates(true);

std::cout << "\nTesting: " << optimizer.name() << " with " << function.name() << " and active updates" << std::endl;
testFunction(optimizer, function, 10, 1000, 1E-10);
}

BOOST_AUTO_TEST_CASE( CMA_Rosenbrock )
{
Rosenbrock function( 3 );
Expand Down
8 changes: 8 additions & 0 deletions include/shark/Algorithms/DirectSearch/CMA.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,13 @@ class CMA : public AbstractSingleObjectiveOptimizer<RealVector >
return m_numEvaluations;
}

bool isUsingActiveUpdates() const {
return m_useActiveUpdates;
}

void useActiveUpsates(bool useActiveUpdates) {
m_useActiveUpdates = useActiveUpdates;
}

protected:
/// \brief The type of individual used for the CMA
Expand All @@ -231,6 +238,7 @@ class CMA : public AbstractSingleObjectiveOptimizer<RealVector >

bool m_userSetMu; /// <The user set a value via setMu, do not overwrite with default
bool m_userSetLambda; /// <The user set a value via setMu, do not overwrite with default
bool m_useActiveUpdates = false; /// <Indicates if active updates should be applied
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please put that in the constructor, everything should be at one place.

double m_initSigma; ///< use provided initial value of sigma<=0 =>algorithm chooses

RecombinationType m_recombinationType; ///< Stores the recombination type.
Expand Down
9 changes: 9 additions & 0 deletions include/shark/Algorithms/DirectSearch/Individual.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,15 @@ class Individual {
}
};

///\brief Reverse ordering relation by the fitness of the individuals(only single objective)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you could use a lambda function for that, we are trying to replace all those small functors with them.

struct ReverseFitnessOrdering
{
bool operator()(Individual const& individual1, Individual const& individual2)
{
return individual1.unpenalizedFitness() < individual2.unpenalizedFitness();
}
};

/// \brief Default constructor that initializes the individual's attributes to default values.
Individual()
: m_rank(0)
Expand Down
100 changes: 83 additions & 17 deletions src/Algorithms/DirectSearch/CMA.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -237,22 +237,31 @@ void CMA::doInit(


//weighting of the k-best individuals
m_weights.resize(m_mu);
m_weights.resize(m_mu * 2);
RealVector negativeWeights(m_mu, 0.);
switch (m_recombinationType) {
case EQUAL:
for (std::size_t i = 0; i < m_mu; i++)
for (std::size_t i = 0; i < m_mu; i++) {
m_weights(i) = 1;
negativeWeights(i) = -1.;
}
break;
case LINEAR:
for (std::size_t i = 0; i < m_mu; i++)
for (std::size_t i = 0; i < m_mu; i++) {
m_weights(i) = (double)(mu-i);
negativeWeights(i) = static_cast<double>(mu - (m_lambda - i - 1.));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you should not need the static_cast because subtracting a double from mu leads to a double, unlike the line before where an unsigned integral type is assigned to double.

}
break;
case SUPERLINEAR:
for (std::size_t i = 0; i < m_mu; i++)
for (std::size_t i = 0; i < m_mu; i++) {
m_weights(i) = ::log(mu + 0.5) - ::log(1. + i); // eq. (45)
negativeWeights(i) = ::log(mu + 0.5) - ::log(1. + (m_lambda - i - 1.));
}
break;
}
m_weights /= sum(m_weights); // eq. (45)
const double weightSum = sum(m_weights);
m_weights /= weightSum; // eq. (45)
negativeWeights /= weightSum; // Normalize the negative weights.
m_muEff = 1. / sum(sqr(m_weights)); // equal to sum(m_weights)^2 / sum(sqr(m_weights))

// Step size control
Expand All @@ -261,9 +270,17 @@ void CMA::doInit(

m_cC = (4. + m_muEff / m_numberOfVariables) / (m_numberOfVariables + 4. + 2 * m_muEff / m_numberOfVariables); // eq. (47)
m_c1 = 2 / (sqr(m_numberOfVariables + 1.3) + m_muEff); // eq. (48)
double alphaMu = 2.;
double rankMuAlpha = 0.3;//but is it really?
m_cMu = std::min(1. - m_c1, alphaMu * ( rankMuAlpha + m_muEff - 2. + 1./m_muEff) / (sqr(m_numberOfVariables + 2) + alphaMu * m_muEff / 2)); // eq. (49)
m_cMu = std::min(1. - m_c1, 2. * (.25 + m_muEff + 1. / m_muEff - 2.) / (std::pow(m_numberOfVariables + 2., 2.) + 2. * m_muEff / 2.)); // eq. (49)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

replace std::pow(a,2.0) by sqr(a)

also this line should be the exact same as before when substitution alphaMu by 2 and changing rankMuAlpha to 0.25


// Normalize the negative weights
const double negativeWeightSum = sum(negativeWeights);
const double negativeMultiplier = 1. + m_c1 / m_cMu;
negativeWeights /= -negativeWeightSum;
negativeWeights *= negativeMultiplier;

for (std::size_t i = 0; i < negativeWeights.size(); ++i) {
m_weights(2 * m_mu - i - 1) = negativeWeights[i];
}

std::size_t pos = std::min_element(initialValues.begin(),initialValues.end())-initialValues.begin();
m_mean = initialSearchPoints[pos];
Expand Down Expand Up @@ -299,13 +316,7 @@ void CMA::updatePopulation( std::vector<IndividualType> const& offspring ) {

// Covariance matrix update
RealMatrix& C = m_mutationDistribution.covarianceMatrix();
RealMatrix Z( m_numberOfVariables, m_numberOfVariables, 0.0); // matric for rank-mu update
for( std::size_t i = 0; i < m_mu; i++ ) {
noalias(Z) += m_weights( i ) * blas::outer_prod(
selectedOffspring[i].searchPoint() - m_mean,
selectedOffspring[i].searchPoint() - m_mean
);
}

double n = static_cast<double>(m_numberOfVariables);
double expectedChi = std::sqrt( n )*(1. - 1./(4.*n) + 1./(21.*n*n));
double hSigLHS = norm_2( m_evolutionPathSigma ) / std::sqrt(1. - pow((1 - m_cSigma), 2.*(m_counter+1)));
Expand All @@ -314,8 +325,63 @@ void CMA::updatePopulation( std::vector<IndividualType> const& offspring ) {
if(hSigLHS < hSigRHS) hSig = 1.;
double deltaHSig = (1.-hSig*hSig) * m_cC * (2. - m_cC);

m_evolutionPathC = (1. - m_cC ) * m_evolutionPathC + hSig * std::sqrt( m_cC * (2. - m_cC) * m_muEff ) * y; // eq. (42)
noalias(C) = (1.-m_c1 - m_cMu) * C + m_c1 * ( blas::outer_prod( m_evolutionPathC, m_evolutionPathC ) + deltaHSig * C) + m_cMu * 1./sqr( m_sigma ) * Z; // eq. (43)
const double c1a = m_c1 * (1. - (1. - (hSig * hSig)) * m_cC * (2. - m_cC));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please move inside the if, if you don't use it in both branches (or later)


if (m_useActiveUpdates) {
// Copy the weights from as they are altered based on the length of the
// the rejected solutions.
RealVector tempWeights(m_numberOfVariables + 1, 0.);
tempWeights(0) = c1a;
for (std::size_t i = 0; i < m_weights.size(); ++i)
{
tempWeights(i + 1) = m_weights(i) * m_cMu;
}

// Add rejected samples to the selected offspring
std::vector< IndividualType > rejectedOffspring(m_mu);
ElitistSelection<IndividualType::ReverseFitnessOrdering > rejectedSelection;
rejectedSelection(offspring.begin(), offspring.end(), rejectedOffspring.begin(), rejectedOffspring.end());
std::reverse(rejectedOffspring.begin(), rejectedOffspring.end());

selectedOffspring.insert(selectedOffspring.end(), rejectedOffspring.begin(), rejectedOffspring.end());

const double weightSum = sum(tempWeights);

const RealMatrix &B = m_mutationDistribution.eigenVectors();
const RealVector &D = sqrt(max(m_mutationDistribution.eigenValues(), 0));

RealMatrix vectors(m_numberOfVariables + 1, m_numberOfVariables);

// Build the matrix for combined rank-1 and rank-mu updates
row(vectors, 0) = m_evolutionPathC * sqrt(m_c1 / (c1a + 1e-23));
for (int k = 0; k < selectedOffspring.size(); ++k)
{
const unsigned int weightIndex = k + 1;
RealVector normalized = (selectedOffspring[k].searchPoint() - m_mean) / m_sigma;

if (tempWeights[weightIndex] < 0.)
{
const double mahalanobisNorm = sqrt(sum(sqr((trans(B) % normalized) / D)));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const double mahalanobisNorm = norm_2((trans(B) % normalized) / D)

or even simpler

const double mahalanobisNorm = norm_2(selectedOffspring[j].chromosome())

tempWeights[weightIndex] *= static_cast<double>(m_numberOfVariables) / sqr(mahalanobisNorm + 1e-9);
}

row(vectors, weightIndex) = normalized;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

noalias(row(vectors,weightIndex)) = normalized;

we don't want to copy!

}

m_evolutionPathC = (1. - m_cC) * m_evolutionPathC + hSig * (std::sqrt(m_cC * (2. - m_cC) * m_muEff) / m_sigma) * (m - m_mean);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is the same as linke 382. (m-mean)/sigma is y. Please try to unify!

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed it should, I did have some weird error when I did it with y instead, but that might origin from somewhere else. I will try to update it!

noalias(C) = C * (1. - weightSum) + trans(to_diagonal(tempWeights) % vectors) % vectors;
} else {
RealMatrix Z(m_numberOfVariables, m_numberOfVariables, 0.0); // matric for rank-mu update
for (std::size_t i = 0; i < m_mu; i++) {
noalias(Z) += m_weights(i) * blas::outer_prod(
selectedOffspring[i].searchPoint() - m_mean,
selectedOffspring[i].searchPoint() - m_mean
);
}

m_evolutionPathC = (1. - m_cC) * m_evolutionPathC + hSig * std::sqrt(m_cC * (2. - m_cC) * m_muEff) * y; // eq. (42)
noalias(C) = (1. - m_c1 - m_cMu) * C + m_c1 * (blas::outer_prod(m_evolutionPathC, m_evolutionPathC) + deltaHSig * C) + m_cMu * 1. / sqr(m_sigma) * Z; // eq. (43)
}

// Step size update
RealVector CInvY = blas::prod( m_mutationDistribution.eigenVectors(), z ); // C^(-1/2)y = Bz
Expand Down