Skip to content

Commit 7165784

Browse files
author
Emmanuel Benazera
committed
improved profile likelihood computation by solving lower dimensional problems + rotating the space for line search, fixes problems on rastrigin + unit tests, ref #48 
1 parent 4d54653 commit 7165784

File tree

4 files changed

+122
-41
lines changed

4 files changed

+122
-41
lines changed

‎src/eo_matrix.h

+9
Original file line numberDiff line numberDiff line change
@@ -61,4 +61,13 @@ inline void removeElement(dVec &vec, unsigned int k)
6161
vec.conservativeResize(vec.size()-1);
6262
}
6363

64+
inline void addElement(dVec &vec, unsigned int k, const double &xk)
65+
{
66+
if (k >= vec.size()+1)
67+
return;
68+
vec.conservativeResize(vec.size()+1);
69+
std::copy(vec.data()+k,vec.data()+vec.size()-1,vec.data()+k+1);
70+
vec[k] = xk;
71+
}
72+
6473
#endif

‎src/errstats.cc

+58-29
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ namespace libcmaes
4747
//debug
4848

4949
pli le(k,samplesize,parameters._dim,parameters._gp.pheno(x),minfvalue,fup,delta);
50-
50+
5151
errstats<TGenoPheno>::profile_likelihood_search(func,parameters,le,cmasol,k,false,samplesize,fup,delta,maxiters,curve); // positive direction
5252
errstats<TGenoPheno>::profile_likelihood_search(func,parameters,le,cmasol,k,true,samplesize,fup,delta,maxiters,curve); // negative direction
5353

@@ -82,15 +82,16 @@ namespace libcmaes
8282
while(true)
8383
{
8484
// get a new xk point.
85-
errstats<TGenoPheno>::take_linear_step(func,parameters,k,minfvalue,fup,delta,n,linit,d,x);
85+
errstats<TGenoPheno>::take_linear_step(func,parameters,k,minfvalue,fup,delta,n,linit,cmasol.eigenvectors(),d,x);
8686
linit = false;
8787

8888
//debug
8989
//std::cout << "new xk point: " << x.transpose() << std::endl;
9090
//debug
9191

9292
// minimize.
93-
CMASolutions ncitsol = errstats<TGenoPheno>::optimize_pk(func,parameters,citsol,k,x[k]);
93+
dVec nx;
94+
CMASolutions ncitsol = errstats<TGenoPheno>::optimize_pk(func,parameters,citsol,k,x,nx);
9495
if (ncitsol._run_status < 0)
9596
{
9697
LOG(ERROR) << "profile likelihood linesearch: optimization error " << ncitsol._run_status << " -- " << ncitsol.status_msg() << std::endl;
@@ -117,13 +118,12 @@ namespace libcmaes
117118
}
118119
else // update current point and solution.
119120
{
120-
//citsol = ncitsol; // XXX: doesn't work well, restarting from min works better.
121-
x = ncitsol.best_candidate().get_x_dvec();
121+
x = nx;
122122
nminfvalue = ncitsol.best_candidate().get_fvalue();
123123
}
124124

125125
// store points.
126-
dVec phenobx = parameters._gp.pheno(ncitsol.best_candidate().get_x_dvec());
126+
dVec phenobx = parameters._gp.pheno(nx);
127127
if (curve)
128128
{
129129
le._fvaluem[samplesize+sign*(1+i)] = ncitsol.best_candidate().get_fvalue();
@@ -148,7 +148,6 @@ namespace libcmaes
148148
}
149149
return;
150150
}
151-
//std::cout << "i=" << i << std::endl;
152151
++i;
153152
if (curve && i == samplesize)
154153
break;
@@ -166,6 +165,7 @@ namespace libcmaes
166165
const double &delta,
167166
const int &n,
168167
const bool &linit,
168+
const dMat &eigenve,
169169
double &d,
170170
dVec &x)
171171
{
@@ -179,18 +179,20 @@ namespace libcmaes
179179
//std::cerr << "d=" << d << " / fdiff=" << fdiff << " / xtmpk=" << xtmp[k] << std::endl;
180180
//debug
181181

182+
dVec dv = dVec::Zero(parameters.dim());
183+
dv[k] = d;
182184
int i = n;
183185
if (fdiff > fup + fdelta) // above
184186
{
185-
//std::cerr << "above\n";
186-
x[k] -= d;
187+
x -= eigenve.transpose() * dv;
187188
while (fdiff > fup + fdelta
188189
&& i > 0
189190
&& phenoxtmp[k] >= parameters._gp.get_boundstrategy().getPhenoLBound(k)
190191
&& phenoxtmp[k] <= parameters._gp.get_boundstrategy().getPhenoUBound(k))
191192
{
192193
d /= 2.0;
193-
xtmp[k] = x[k] + d;
194+
dv[k] = d;
195+
xtmp = x + eigenve.transpose() * dv; // search in rotated space
194196
phenoxtmp = parameters._gp.pheno(xtmp);
195197
fvalue = func(phenoxtmp.data(),xtmp.size());
196198
fdiff = fabs(fvalue - minfvalue);
@@ -205,14 +207,14 @@ namespace libcmaes
205207
}
206208
else // below
207209
{
208-
//std::cerr << "below\n";
209210
while (fdiff < fup - fdelta
210211
&& i > 0
211212
&& phenoxtmp[k] >= parameters._gp.get_boundstrategy().getPhenoLBound(k)
212213
&& phenoxtmp[k] <= parameters._gp.get_boundstrategy().getPhenoUBound(k))
213214
{
214215
d *= 2.0;
215-
xtmp[k] = x[k] + d;
216+
dv[k] = d;
217+
xtmp = x + eigenve.transpose() * dv; // search in rotated space
216218
phenoxtmp = parameters._gp.pheno(xtmp);
217219
fvalue = func(phenoxtmp.data(),xtmp.size());
218220
fdiff = fabs(fvalue - minfvalue);
@@ -234,30 +236,56 @@ namespace libcmaes
234236
const CMAParameters<TGenoPheno> &parameters,
235237
const CMASolutions &cmasol,
236238
const std::vector<int> &k,
237-
const std::vector<double> &vk)
239+
const std::vector<double> &vk,
240+
const dVec &x0,
241+
dVec &nx)
238242
{
239-
CMASolutions ncmasol = cmasol;
240-
CMAParameters<TGenoPheno> nparameters = parameters;
241-
nparameters._quiet = true;
242-
ncmasol._sigma = nparameters._sigma_init;
243+
dVec rx0 = x0;
244+
double sigma = 0.0;
243245
for (size_t i=0;i<k.size();i++)
244246
{
245-
nparameters.set_fixed_p(k[i],vk[i]);
246-
nparameters._sigma_init = ncmasol._sigma = std::max(ncmasol._sigma,fabs(cmasol.best_candidate().get_x_dvec_ref()[k[i]]-vk[i])); // XXX: here sigma as max distance to best candidate in parameter space. There might be better selection schemes.
247+
sigma = std::max(sigma,fabs(cmasol.best_candidate().get_x_dvec_ref()[k[i]]-vk[i]));
248+
removeElement(rx0,k[i]);
247249
}
248-
return cmaes<TGenoPheno>(func,nparameters,CMAStrategy<CovarianceUpdate,TGenoPheno>::_defaultPFunc,nullptr,ncmasol); // XXX: we're not explicitely setting the initial covariance because current set of experiments shows that it doesn't work well.
250+
CMAParameters<TGenoPheno> nparameters(rx0.size(),rx0.data(),sigma,parameters.lambda());
251+
nparameters.set_ftarget(parameters.get_ftarget());
252+
253+
FitFunc rfunc = [func,k,vk](const double *x, const int N)
254+
{
255+
std::vector<double> nx;
256+
nx.reserve(N+1);
257+
for (size_t j=0;j<k.size();j++)
258+
{
259+
for (int i=0;i<N+1;i++)
260+
{
261+
if (i < k[j])
262+
nx.push_back(x[i]);
263+
else if (i == k[j])
264+
nx.push_back(vk[j]);
265+
else nx.push_back(x[i-1]);
266+
}
267+
}
268+
return func(&nx.front(),N+1);
269+
};
270+
271+
CMASolutions cms = cmaes<TGenoPheno>(rfunc,nparameters);
272+
nx = cms.best_candidate().get_x_dvec();
273+
for (size_t i=0;i<k.size();i++)
274+
addElement(nx,k[i],vk[i]);
275+
return cms;
249276
}
250277

251278
template <class TGenoPheno>
252279
CMASolutions errstats<TGenoPheno>::optimize_pk(FitFunc &func,
253280
const CMAParameters<TGenoPheno> &parameters,
254281
const CMASolutions &cmasol,
255282
const int &k,
256-
const double &vk)
283+
const dVec &vk,
284+
dVec &nx)
257285
{
258286
std::vector<int> tk = {k};
259-
std::vector<double> tvk = {vk};
260-
return errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,tk,tvk);
287+
std::vector<double> tvk = {vk[k]};
288+
return errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,tk,tvk,vk,nx);
261289
}
262290

263291
template <class TGenoPheno>
@@ -423,13 +451,14 @@ namespace libcmaes
423451
double tlf = ftol*fup;
424452
int nfcn = 0;
425453
unsigned int maxitr = 15, ipt = 0;
454+
dVec nx;
426455

427456
//debug
428457
//std::cout << "aminsv=" << aminsv << " / aim=" << aim << " / tla=" << tla << " / tlf=" << tlf << std::endl;
429458
//debug
430459

431-
// get a first optimized point.
432-
CMASolutions cmasol1 = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmid);
460+
// get a first optimized point.
461+
CMASolutions cmasol1 = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmid,parameters.get_x0min(),nx);
433462
alsb[0] = 0.0;
434463
flsb[0] = cmasol1.best_candidate().get_fvalue();
435464
flsb[0] = std::max(flsb[0],aminsv+0.1*fup);
@@ -448,7 +477,7 @@ namespace libcmaes
448477
else if (aopt < -0.5)
449478
aopt = 0.5;
450479
std::vector<double> pmiddir2 = {pmid[0]+aopt*pdir[0],pmid[1]+aopt*pdir[1]};
451-
CMASolutions cmasol2 = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmiddir2);
480+
CMASolutions cmasol2 = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmiddir2,parameters.get_x0min(),nx);
452481
alsb[1] = aopt;
453482
flsb[1] = cmasol2.best_candidate().get_fvalue();
454483
nfcn += cmasol2._nevals;
@@ -477,7 +506,7 @@ namespace libcmaes
477506
flsb[0] = flsb[1];
478507
aopt = alsb[0] + 0.2*(it+1);
479508
std::vector<double> pmidt = {pmid[0]+aopt*pdir[0],pmid[1]+aopt*pdir[1]};
480-
CMASolutions cmasolt = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmidt);
509+
CMASolutions cmasolt = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmidt,parameters.get_x0min(),nx);
481510
alsb[1] = aopt;
482511
flsb[1] = cmasolt.best_candidate().get_fvalue();
483512
dfda = (flsb[1]-flsb[0])/(alsb[1]-alsb[0]);
@@ -517,7 +546,7 @@ namespace libcmaes
517546

518547
// get third point.
519548
std::vector<double> pmiddir3 = {pmid[0]+aopt*pdir[0],pmid[1]+aopt*pdir[1]};
520-
CMASolutions cmasol3 = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmiddir3);
549+
CMASolutions cmasol3 = errstats<TGenoPheno>::optimize_vpk(func,parameters,cmasol,par,pmiddir3,parameters.get_x0min(),nx);
521550
alsb[2] = aopt;
522551
flsb[2] = cmasol3.best_candidate().get_fvalue();
523552
nfcn += cmasol3._nevals;
@@ -609,7 +638,7 @@ namespace libcmaes
609638
//debug
610639

611640
std::vector<double> vxk = {xstart(par[0]),xstart(par[1])};
612-
CMASolutions ncitsol = errstats<TGenoPheno>::optimize_vpk(func,parameters,citsol,par,vxk);
641+
CMASolutions ncitsol = errstats<TGenoPheno>::optimize_vpk(func,parameters,citsol,par,vxk,parameters.get_x0min(),nx);
613642
if (ncitsol._run_status < 0)
614643
{
615644
LOG(WARNING) << "contour linesearch: optimization error " << ncitsol._run_status << std::endl;

‎src/errstats.h

+9-2
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,9 @@ namespace libcmaes
9696
* @param fup the function deviation for which to compute the profile likelihood
9797
* @param delta tolerance around fvalue + fup for which to compute the profile likelihood
9898
* @param linit whether this is the first linesearch call
99+
* @param eigenve eigenvectors
100+
* @param d step
101+
* @param x vector on the line
99102
*/
100103
static void take_linear_step(FitFunc &func,
101104
const CMAParameters<TGenoPheno> &parameters,
@@ -105,6 +108,7 @@ namespace libcmaes
105108
const double &delta,
106109
const int &n,
107110
const bool &linit,
111+
const dMat &eigenve,
108112
double &d,
109113
dVec &x);
110114

@@ -122,7 +126,9 @@ namespace libcmaes
122126
const CMAParameters<TGenoPheno> &parameters,
123127
const CMASolutions &cmasol,
124128
const std::vector<int> &k,
125-
const std::vector<double> &vk);
129+
const std::vector<double> &vk,
130+
const dVec &x0,
131+
dVec &nx);
126132

127133
/**
128134
* \brief optimizes an objective function while fixing the value of parameters in dimension k
@@ -137,7 +143,8 @@ namespace libcmaes
137143
const CMAParameters<TGenoPheno> &parameters,
138144
const CMASolutions &cmasol,
139145
const int &k,
140-
const double &vk);
146+
const dVec &vk,
147+
dVec &nx);
141148

142149
/*- contour -*/
143150
public:

‎tests/ut-errstats.cc

+46-10
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,44 @@
2525

2626
using namespace libcmaes;
2727

28+
TEST(eomatrix,removeElement)
29+
{
30+
for (int k=0;k<10;k++)
31+
{
32+
dVec x = dVec::Random(10);
33+
dVec xp = x;
34+
removeElement(xp,k);
35+
ASSERT_EQ(9,xp.size());
36+
for (int i=0;i<10;i++)
37+
{
38+
if (i<k)
39+
ASSERT_EQ(x[i],xp[i]);
40+
else if (i>k)
41+
ASSERT_EQ(x[i],xp[i-1]);
42+
}
43+
}
44+
}
45+
46+
TEST(eomatrix,addElement)
47+
{
48+
for (int k=0;k<11;k++)
49+
{
50+
dVec x = dVec::Random(10);
51+
dVec xp = x;
52+
addElement(xp,k,2.5);
53+
ASSERT_EQ(11,xp.size());
54+
for (int i=0;i<1;i++)
55+
{
56+
if (i<k)
57+
ASSERT_EQ(x[i],xp[i]);
58+
else if (i == k)
59+
ASSERT_EQ(2.5,xp[i]);
60+
else if (i>k)
61+
ASSERT_EQ(x[i+1],xp[i]);
62+
}
63+
}
64+
}
65+
2866
TEST(rearrangecmasol,reset_as_fixed)
2967
{
3068
FitFunc fsphere = [](const double *x, const int N)
@@ -45,8 +83,6 @@ TEST(rearrangecmasol,reset_as_fixed)
4583
ASSERT_EQ(9,copsols.cov().rows());
4684
ASSERT_EQ(9,copsols.cov().cols());
4785
ASSERT_EQ(9,copsols.xmean().size());
48-
/*std::cout << cmasols.xmean().transpose() << std::endl;
49-
std::cout << copsols.xmean().transpose() << std::endl;*/
5086
CMAParameters<> copparams = cmaparams;
5187
cmaparams.reset_as_fixed(6);
5288

@@ -67,13 +103,13 @@ TEST(optimize,optimize_fixed_p)
67103
CMAParameters<> cmaparams(dim,&x0.front(),sigma);
68104
cmaparams.set_quiet(true);
69105
CMASolutions cmasols = cmaes<>(fsphere,cmaparams);
70-
/*cmaparams.set_fixed_p(6,0.1);
71-
CMASolutions cmaksols = cmaes<>(fsphere,cmaparams);*/
72-
CMASolutions cmaksols = errstats<>::optimize_pk(fsphere,cmaparams,cmasols,6,0.1);
106+
dVec nx;
107+
CMASolutions cmaksols = errstats<>::optimize_pk(fsphere,cmaparams,cmasols,6,cmasols.xmean(),nx);
73108
std::cout << "iter: " << cmaksols.niter() << std::endl;
74109
std::cout << "run status: " << cmaksols.run_status() << std::endl;
75110
ASSERT_EQ(TOLHISTFUN,cmaksols.run_status());
76-
ASSERT_EQ(10,cmaksols.best_candidate().get_x_dvec().size());
111+
ASSERT_EQ(9,cmaksols.best_candidate().get_x_dvec().size());
112+
ASSERT_EQ(10,nx.size());
77113
std::cout << "fvalue: " << cmaksols.best_candidate().get_fvalue() << std::endl;
78114
std::cout << "x: " << cmaksols.best_candidate().get_x_dvec().transpose() << std::endl;
79115
}
@@ -100,8 +136,8 @@ TEST(pl,profile_likelihood_nocurve)
100136
pli le = errstats<>::profile_likelihood(fsphere,cmaparams,cmasols,k,false,samplesize,fup);
101137
std::cout << "le fvalue: " << le.get_fvaluem().transpose() << std::endl;
102138
std::cout << "le xm: " << le.get_xm() << std::endl;
103-
EXPECT_FLOAT_EQ(-0.31640971,le.get_min());
104-
EXPECT_FLOAT_EQ(0.31640971,le.get_max());
139+
EXPECT_FLOAT_EQ(-0.32090676,le.get_min());
140+
EXPECT_FLOAT_EQ(0.32090676,le.get_max());
105141
}
106142

107143
TEST(pl,profile_likelihood_curve)
@@ -128,6 +164,6 @@ TEST(pl,profile_likelihood_curve)
128164
std::cout << "le xm: " << le.get_xm() << std::endl;
129165
int mini, maxi;
130166
std::pair<double,double> mm = le.getMinMax(0.1,mini,maxi);
131-
EXPECT_FLOAT_EQ(-0.311827,mm.first);
132-
EXPECT_FLOAT_EQ(0.311827,mm.second);
167+
EXPECT_FLOAT_EQ(-0.30449873,mm.first);
168+
EXPECT_FLOAT_EQ(0.30449873,mm.second);
133169
}

0 commit comments

Comments
 (0)