From 4c37ffdf1215c5b66a529453aa608bec72c3c7db Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Mon, 26 Oct 2026 12:52:04 +0100 Subject: [PATCH] optimize: added Powell and Nelder-Mead method --- optimize.cc | 302 ++++++++++++++++++++++++++++++++++++++++++++++++++++ optimize.h | 50 +++++++++ t.cc | 16 ++- 3 files changed, 363 insertions(+), 5 deletions(-) diff --git a/optimize.cc b/optimize.cc index 1532c2a..3dcd974 100644 --- a/optimize.cc +++ b/optimize.cc @@ -233,6 +233,112 @@ return optresult(f,converged); } +//function along the line x0 + t*u for linemin_brent, keeps the best point and counts the evaluations +template +class linefunction { +public: +optevaluator &e; +const NRVec &x0; +const NRVec &u; +int maxiter; //evaluation budget +int nev; //evaluations done +T tbest; +T fbest; +linefunction(optevaluator &ee, const NRVec &xx0, const T f0, const NRVec &uu, const int maxit) : e(ee), x0(xx0), u(uu), maxiter(std::min(maxit,ee.remaining())), nev(0), tbest(0), fbest(f0) {}; +NRVec point(const T t) const {NRVec x(x0); x.axpy(t,u); return x;} +T operator()(const T t) {T f=e(point(t)); ++nev; if(f=maxiter;} +}; + + +template +T linemin_brent(optevaluator &e, NRVec &x, T &f, const NRVec &u, const T t1, const T tol, const T abstol, const int maxiter, const int verbose) +{ +const T golden=1.618033988749895; //expansion factor of the bracketing +const T cgold=0.3819660112501051; //golden section fraction (3-sqrt(5))/2 +linefunction phi(e,x,f,u,maxiter); + +//bracketing: a,b,c with f(b)<=f(a), f(b)<=f(c) +T a=0, fa=f; +T b=t1, fb; +T c=0, fc=0; +bool bracketed=false; +if(phi.exhausted()) return 0; +fb=phi(b); +if(fb>fa) //uphill, try the opposite direction + { + if(phi.exhausted()) return 0; + c= -t1; fc=phi(c); + if(fc>=fa) {std::swap(a,b); std::swap(fa,fb); bracketed=true;} //minimum between -t1 and t1, middle point is 0 + else {b=c; fb=fc;} + } +if(!bracketed) //expand downhill from a through b + { + for(;;) + { + if(phi.exhausted()) break; + c= b+golden*(b-a); fc=phi(c); + if(fc>=fb) {bracketed=true; break;} + a=b; fa=fb; b=c; fb=fc; + } + } +if(verbose>1) std::cout <<"linemin bracket "<::epsilon()*std::abs(t1)); //nonzero to guarantee termination + while(!phi.exhausted()) + { + T mid=(lo+hi)/2; + T tol1=tol*std::abs(xm)+tolabs; + T tol2=2*tol1; + if(std::abs(xm-mid) <= tol2-(hi-lo)/2) break; //converged + bool goldenstep=true; + if(std::abs(dprev)>tol1) //try parabolic interpolation through xm, w, v + { + T r=(xm-w)*(fx-fv); + T q=(xm-v)*(fx-fw); + T p=(xm-v)*q-(xm-w)*r; + q=2*(q-r); + if(q>0) p= -p; else q= -q; + if(std::abs(p)q*(lo-xm) && pxm ? tol1 : -tol1; + goldenstep=false; + } + } + if(goldenstep) {dprev= xm>=mid ? lo-xm : hi-xm; d=cgold*dprev;} + T t= std::abs(d)>=tol1 ? xm+d : xm+(d>0 ? tol1 : -tol1); + T ft=phi(t); + if(ft<=fx) + { + if(t>=xm) lo=xm; else hi=xm; + v=w; fv=fw; w=xm; fw=fx; xm=t; fx=ft; + } + else + { + if(t1) std::cout <<"linemin t= "< optresult linesearch_minimize(const char *name, linesearch_direction &dir, const OptProblem *o, NRVec &x, const linesearch_hyperparameters &hp, hyperparameters *achieved) { @@ -442,6 +548,199 @@ return linesearch_minimize(names[hp.my_beta],dir,o,x,hp,achieved); +//simplex of n+1 vertices with their function values +template +class NM_simplex { +public: +int n; +NRVec > v; //vertices +NRVec f; //function values +NRVec sum; //sum of the vertices +int best, second, worst; //indices of the lowest, second highest and highest function value +NM_simplex(const NRVec &x0, const T h, optevaluator &e); +void order(); +void sumup() {sum=v[0]; for(int i=1; i<=n; ++i) sum+=v[i];} +void replace(const int i, const NRVec &x, const T fx) {sum-=v[i]; sum+=x; v[i]=x; f[i]=fx;} +NRVec centroid() const {NRVec c(sum); c-=v[worst]; c/=(T)n; return c;} //of all vertices except the worst one +NRVec point(const NRVec &c, const T t) const {NRVec x = c*(1+t); x.axpy(-t,v[worst]); return x;} //c + t*(c - x_worst) +T diameter() const; +void shrink(const T sigma, optevaluator &e); +}; + +//initial simplex x0, x0 + h e_i +template +NM_simplex::NM_simplex(const NRVec &x0, const T h, optevaluator &e) : n(x0.size()), v(n+1), f(n+1), best(0), second(0), worst(0) +{ +v[0]=x0; +f[0]=e(x0); +for(int i=1; i<=n; ++i) + { + v[i]=x0; + v[i].copyonwrite(); + v[i][i-1] += h; + f[i]=e(v[i]); + } +sumup(); +} + +template +void NM_simplex::order() +{ +best=worst=0; +for(int i=1; i<=n; ++i) + { + if(f[i]=f[worst]) worst=i; + } +second=best; +for(int i=0; i<=n; ++i) if(i!=worst && f[i]>=f[second]) second=i; +} + +template +T NM_simplex::diameter() const +{ +T d=0; +for(int i=0; i<=n; ++i) if(i!=best) + { + NRVec dx = v[i]-v[best]; + d=std::max(d,dx.norm()); + } +return d; +} + +//all vertices move towards the best one, as far as the evaluation budget allows +template +void NM_simplex::shrink(const T sigma, optevaluator &e) +{ +for(int i=0; i<=n; ++i) if(i!=best) + { + if(e.exhausted()) break; + NRVec x = v[best]*(1-sigma); + x.axpy(sigma,v[i]); + v[i]=x; + f[i]=e(x); + } +sumup(); +} + + +template +optresult NelderMead_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +const char *name="Nelder-Mead"; +NelderMead_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("Nelder-Mead cannot search for saddle points"); +const int verbose=hp.my_verbose; +const int n=x.size(); +T alpha=hp.my_alpha, gamma=hp.my_gamma, rho=hp.my_rho, sigma=hp.my_sigma; +if(hp.my_adaptive) {alpha=1; gamma=1+(T)2/n; rho=0.75-(T)0.5/n; sigma=1-(T)1/n;} +optevaluator e(o,hp.my_maxfunc); +optstate st(name,hp,e); +NM_simplex s(x,hp.my_initstep,e); +const T huge=std::numeric_limits::max(); + +for(;; ++st.iter) + { + s.order(); + T spread=s.f[s.worst]-s.f[s.best]; + st.step(spread, spread c=s.centroid(); + NRVec xr=s.point(c,alpha); + T fr=e(xr); + if(fr xe=s.point(c,alpha*gamma); + T fe=e(xe); + if(fe xc=s.point(c, outside ? alpha*rho : -rho); + T fc=e(xc); + if(outside ? fc<=fr : fc1) std::cout < +optresult Powell_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +const char *name="Powell"; +Powell_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("Powell cannot search for saddle points"); +const int verbose=hp.my_verbose; +const T abstol=0.01*hp.my_x_thr; //line minima have to be more precise than the requested accuracy in x +const int n=x.size(); +optevaluator e(o,hp.my_maxfunc); +optstate st(name,hp,e); +T f=e(x); + +NRVec > u(n); //normalized directions +for(int i=0; i((T)0,n); u[i][i]=1;} +NRVec tstep(hp.my_initstep,n); //trial steps of the line minimizations, from the last step along each direction + +for(;; ++st.iter) + { + if(st.stop(f)) break; + NRVec x0(x); + T f0=f; + T delta=0; //largest decrease along a single direction + int ibig=0; + for(int i=0; idelta) {delta=fprev-f; ibig=i;} + } + + //new direction x-x0 replaces the one of the largest decrease, if Powell's test indicates it is worth + NRVec unew = x-x0; + T unorm=unew.norm(); + if(unorm>0 && !e.exhausted()) + { + NRVec xe = x+unew; //extrapolated point + T fe=e(xe); + T a=f0-f-delta; + T b=f0-fe; + if(fe1) std::cout < dx = x-x0; + st.step(f0-f,dx.norm()); + } + +return st.finish(achieved); +} + + + //force instantiation template double linesearch_wolfe(optevaluator &e, NRVec &x, double &f, NRVec &g, const NRVec &p, const double alpha1, const linesearch_hyperparameters &lp); template class optstate; @@ -449,6 +748,9 @@ template class testproblem; template class BFGS_Optimizer; template class LBFGS_Optimizer; template class CG_Optimizer; +template class NelderMead_Optimizer; +template class Powell_Optimizer; +template double linemin_brent(optevaluator &e, NRVec &x, double &f, const NRVec &u, const double t1, const double tol, const double abstol, const int maxiter, const int verbose); }//namespace diff --git a/optimize.h b/optimize.h index 76d57fe..6c854dc 100644 --- a/optimize.h +++ b/optimize.h @@ -172,6 +172,13 @@ template T linesearch_wolfe(optevaluator &e, NRVec &x, T &f, NRVec &g, const NRVec &p, const T alpha1, const linesearch_hyperparameters &lp); +//derivative-free line minimization of f(x + t*u): bracketing and Brent's method (Brent, Algorithms for Minimization without Derivatives, 1973, chap. 5) +//on input x,f at the starting point, on output x,f at the minimum found; returns t, 0 if no decrease was found (x,f then unchanged) +//t1 = initial trial step, tol and abstol = relative and absolute tolerance on t, maxiter = max. no. of function evaluations +template +T linemin_brent(optevaluator &e, NRVec &x, T &f, const NRVec &u, const T t1, const T tol, const T abstol, const int maxiter, const int verbose=0); + + //search direction of a line-search method, used by the common driver linesearch_minimize() template class linesearch_direction { @@ -261,6 +268,49 @@ virtual ~CG_Optimizer() {}; }; +//Nelder-Mead downhill simplex (Lagarias et al., SIAM J. Optim. 9, 112 (1998)), derivative-free +//convergence: f_thr is used for the spread of function values in the simplex, x_thr for its diameter, both have to be satisfied +template +class NelderMead_hyperparameters : public hyperparameters { +public: +T my_initstep; //edge length of the initial simplex x0 + initstep*e_i +bool my_adaptive; //dimension-dependent coefficients (Gao & Han, Comput. Optim. Appl. 51, 259 (2012)), otherwise the following ones +T my_alpha; //reflection +T my_gamma; //expansion +T my_rho; //contraction +T my_sigma; //shrink +NelderMead_hyperparameters() : my_initstep(0.1), my_adaptive(true), my_alpha(1), my_gamma(2), my_rho(0.5), my_sigma(0.5) {this->my_maxstep=5000; this->my_f_thr=1e-8; this->my_x_thr=1e-6;}; +virtual void setfrom(const hyperparameters &p) {const NelderMead_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} +}; + +template +class NelderMead_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~NelderMead_Optimizer() {}; +}; + + +//Powell's conjugate direction method (Powell, Comput. J. 7, 155 (1964)), derivative-free +//convergence: f_thr and x_thr are used for the change of f and x in one sweep over all directions, both have to be satisfied +template +class Powell_hyperparameters : public hyperparameters { +public: +T my_initstep; //initial trial step of the line minimizations +T my_linmin_tol; //relative tolerance of the line minimizations, the absolute one is 0.01*x_thr +int my_maxlinemin; //max. no. of function evaluations in one line minimization +Powell_hyperparameters() : my_initstep(0.1), my_linmin_tol(1e-4), my_maxlinemin(50) {this->my_maxstep=1000; this->my_f_thr=1e-10; this->my_x_thr=1e-6;}; +virtual void setfrom(const hyperparameters &p) {const Powell_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} +}; + +template +class Powell_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~Powell_Optimizer() {}; +}; + + }//namespace #endif /* _OPTIMIZE_H_ */ diff --git a/t.cc b/t.cc index ce55aa1..b3b5f83 100644 --- a/t.cc +++ b/t.cc @@ -4845,22 +4845,28 @@ x.randomize(1.); testproblem problem; //BFGS_hyperparameters param; //tested //LBFGS_hyperparameters param; //tested -CG_hyperparameters param; +//CG_hyperparameters param; //tested +//NelderMead_hyperparameters param; //tested +Powell_hyperparameters param; param.my_maxstep=100; param.my_f_thr=1e-6; param.my_x_thr=1e-3; -param.my_g_thr=1e-3; +//param.my_g_thr=1e-3; //gradient methods //param.my_update=QN_BFGS; //tested //param.my_update=QN_DFP; //tested //param.my_update=QN_SR1; //tested cout << "initial " < o; //tested //LBFGS_Optimizer o; //tested -CG_Optimizer o; -gradient_hyperparameters achieved; +//CG_Optimizer o; //tested +//NelderMead_Optimizer o; //tested +Powell_Optimizer o; +//gradient_hyperparameters achieved; //gradient methods +hyperparameters achieved; optresult res=o.minimize(&problem,x,¶m,&achieved); cout << "optimized " <