diff --git a/optimize.cc b/optimize.cc index ebb5a6b..c3dd1b5 100644 --- a/optimize.cc +++ b/optimize.cc @@ -1124,7 +1124,7 @@ NRSMat A=J.transposedtimes(); //J^T J NRVec g=r*J; //J^T r, the gradient of f NRVec diag(n); for(int i=0; i &)A)(i,i); -T mu=hp.my_tau*diag.amax(), nu=2; +T mu=hp.my_scaling ? hp.my_tau : hp.my_tau*diag.amax(), nu=2; for(; budget; ++st.iter) { @@ -1163,7 +1163,12 @@ for(; budget; ++st.iter) { mu*=nu; nu*=2; const T eps=std::numeric_limits::epsilon(); - if(h.norm()<=10*eps*(x.norm()+eps)) //no further improvement at machine precision (cf. Madsen, Nielsen & Tingleff eq. 3.16) + //no further improvement at machine precision (cf. Madsen, Nielsen & Tingleff eq. 3.16), tested per component to be scale invariant: + //h_i^2 <= (10 eps)^2 (x_i^2+eps^2) for all i (with |h|<=10 eps |x| badly scaled parameters stopped too early in float) + NRVec hh(h); hh*=h; + NRVec xx(x); xx*=x; xx+=eps*eps; + hh/=xx; + if(hh.amax()<=100*eps*eps) { if(verbose) std::cout < firstorder_minimize(const char *name, firstorder_step & template optresult trustregion_minimize(const char *name, trustregion_step &stepper, const OptProblem *o, NRVec &x, const trustregion_hyperparameters &hp, hyperparameters *achieved); \ INSTANTIZE(double) +INSTANTIZE(float) }//namespace diff --git a/optimize.h b/optimize.h index 191939c..d4269cf 100644 --- a/optimize.h +++ b/optimize.h @@ -34,6 +34,37 @@ namespace LA { +//precision-dependent default steps and thresholds of the optimizers; only double and float are supported (the double values are those used so far) +template struct opt_precision; +template<> struct opt_precision + { + static double fd_central_step() {return 1e-6;} //central differences (FiniteDifference, ParameterShift), optimum ~ eps^(1/3) + static double fd_forward_step() {return 1e-7;} //forward differences (OptLeastSquares, LM), optimum ~ eps^(1/2) + static double g_thr() {return 1e-6;} //gradient convergence + static double sy_thr() {return 1e-12;} //curvature condition of the quasi-Newton updates + static double sr1_thr() {return 1e-8;} //SR1 skipping rule + static double update_thr() {return 1e-8;} //trust-region Hessian update skipping rule + static double nm_f_thr() {return 1e-8;} //Nelder-Mead stopping + static double nm_x_thr() {return 1e-6;} + static double powell_f_thr() {return 1e-10;} //Powell stopping + static double powell_x_thr() {return 1e-6;} + static double rhoend() {return 1e-6;} //final COBYLA trust radius + }; +template<> struct opt_precision + { + static float fd_central_step() {return 5e-3f;} + static float fd_forward_step() {return 3e-4f;} + static float g_thr() {return 1e-3f;} + static float sy_thr() {return 1e-6f;} + static float sr1_thr() {return 1e-4f;} + static float update_thr() {return 1e-4f;} + static float nm_f_thr() {return 1e-5f;} + static float nm_x_thr() {return 1e-4f;} + static float powell_f_thr() {return 1e-6f;} + static float powell_x_thr() {return 1e-4f;} + static float rhoend() {return 1e-4f;} + }; + //base classes with purely virtual functions for multiple inheritance use //optimization problem - provides function and optional gradient, in multiply inherited classess will store auxiliary data for the function @@ -102,7 +133,7 @@ template class OptLeastSquares : public OptProblem, public OptResiduals { public: T fd_step; //relative step of the forward-difference Jacobian -OptLeastSquares(const T h=1e-7) : fd_step(h) {}; +OptLeastSquares(const T h=opt_precision::fd_forward_step()) : fd_step(h) {}; void jacobian_fd(const NRVec &x, const NRVec &r, NRMat &J) const; //Jacobian by forward differences, r = residuals at x virtual T funcgrad(const NRVec &x, NRVec *grad) const; virtual ~OptLeastSquares() {}; @@ -164,7 +195,7 @@ class ParameterShift : public OptProblemAdapter { public: const OptTrigonometric *trig; T h; -ParameterShift(const OptProblem *pp, const T hh=1e-6); +ParameterShift(const OptProblem *pp, const T hh=opt_precision::fd_central_step()); virtual T funcgrad(const NRVec &x, NRVec *grad) const; }; @@ -173,7 +204,7 @@ template class FiniteDifference : public OptProblemAdapter { public: T h; -FiniteDifference(const OptProblem *pp, const T hh=1e-6) : OptProblemAdapter(pp), h(hh) {}; +FiniteDifference(const OptProblem *pp, const T hh=opt_precision::fd_central_step()) : OptProblemAdapter(pp), h(hh) {}; virtual T funcgrad(const NRVec &x, NRVec *grad) const; }; @@ -207,7 +238,7 @@ template class gradient_hyperparameters : public hyperparameters { public: T my_g_thr; //gradient norm convergence threshold -gradient_hyperparameters() : my_g_thr(1e-6) {}; +gradient_hyperparameters() : my_g_thr(opt_precision::g_thr()) {}; virtual void setfrom(const hyperparameters &p) {const gradient_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} virtual bool converged(const optstate &s) const {return (s.gnorm>=0 && s.gnorm::converged(s);} virtual void setreached(const optstate &s) {hyperparameters::setreached(s); my_g_thr=s.gnorm;} @@ -379,7 +410,7 @@ template class quasinewton_hyperparameters : public linesearch_hyperparameters { public: T my_sy_thr; //skip the BFGS-type update if s.y <= sy_thr*|s|*|y| (curvature condition violated numerically) -quasinewton_hyperparameters() : my_sy_thr(1e-12) {}; +quasinewton_hyperparameters() : my_sy_thr(opt_precision::sy_thr()) {}; virtual void setfrom(const hyperparameters &p) {const quasinewton_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else linesearch_hyperparameters::setfrom(p);} }; @@ -392,7 +423,7 @@ class BFGS_hyperparameters : public quasinewton_hyperparameters { public: QN_update my_update; //BFGS, DFP or SR1 (Nocedal & Wright eqs. 6.17, 6.15, 6.25); DFP needs an accurate line search (c2~0.1), otherwise it can stall T my_sr1_thr; //SR1: skip the update if |(s-Hy).y| <= sr1_thr*|s-Hy|*|y| (cf. Nocedal & Wright eq. 6.26) -BFGS_hyperparameters() : my_update(QN_BFGS), my_sr1_thr(1e-8) {}; +BFGS_hyperparameters() : my_update(QN_BFGS), my_sr1_thr(opt_precision::sr1_thr()) {}; virtual void setfrom(const hyperparameters &p) {const BFGS_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else quasinewton_hyperparameters::setfrom(p);} }; @@ -485,7 +516,7 @@ T my_eta; //a step of a minimization is accepted if (actual reduction)/(predicte TR_hessian my_hessian; int my_hessian_recompute; //exact Hessian computed every k-th iteration and updated in between; 1 = always exact T my_update_thr; //threshold for skipping the Hessian updates (curvature condition of BFGS, denominator of SR1) -trustregion_hyperparameters() : my_delta0(1), my_deltamax(1e3), my_eta(1e-4), my_hessian(TR_AUTO), my_hessian_recompute(1), my_update_thr(1e-8) {}; +trustregion_hyperparameters() : my_delta0(1), my_deltamax(1e3), my_eta(1e-4), my_hessian(TR_AUTO), my_hessian_recompute(1), my_update_thr(opt_precision::update_thr()) {}; virtual void setfrom(const hyperparameters &p) {const trustregion_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else gradient_hyperparameters::setfrom(p);} }; @@ -553,7 +584,7 @@ public: T my_tau; //initial damping mu = tau*max diag(J^T J) bool my_scaling; //Marquardt scaling D = diag(J^T J) T my_fd_step; //relative step of the forward-difference Jacobian -LM_hyperparameters() : my_tau(1e-3), my_scaling(false), my_fd_step(1e-7) {}; +LM_hyperparameters() : my_tau(1e-3), my_scaling(false), my_fd_step(opt_precision::fd_forward_step()) {}; virtual void setfrom(const hyperparameters &p) {const LM_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else gradient_hyperparameters::setfrom(p);} }; @@ -576,7 +607,7 @@ 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;}; +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=opt_precision::nm_f_thr(); this->my_x_thr=opt_precision::nm_x_thr();}; virtual void setfrom(const hyperparameters &p) {const NelderMead_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} }; @@ -596,7 +627,7 @@ 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;}; +Powell_hyperparameters() : my_initstep(0.1), my_linmin_tol(1e-4), my_maxlinemin(50) {this->my_maxstep=1000; this->my_f_thr=opt_precision::powell_f_thr(); this->my_x_thr=opt_precision::powell_x_thr();}; virtual void setfrom(const hyperparameters &p) {const Powell_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} }; @@ -620,7 +651,7 @@ T my_rhobeg; //initial trust region radius, about one tenth of the expected chan T my_rhoend; //final trust region radius, the accuracy of x bool my_flexible_radius; //trust region radius delta >= rho revised by the reduction ratio and rho reduced by the rules of NEWUOA //(Powell, Report DAMTP 2004/NA08, eqs. (7.1), (7.6), (6.16)), so that delta can also increase; false = Powell's original COBYLA with delta = rho -Cobyla_hyperparameters() : my_rhobeg(0.5), my_rhoend(1e-6), my_flexible_radius(true) {this->my_maxstep=100000;}; +Cobyla_hyperparameters() : my_rhobeg(0.5), my_rhoend(opt_precision::rhoend()), my_flexible_radius(true) {this->my_maxstep=100000;}; virtual void setfrom(const hyperparameters &p) {const Cobyla_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} }; diff --git a/testoptimize.cc b/testoptimize.cc index 5e14d60..eee9cf3 100644 --- a/testoptimize.cc +++ b/testoptimize.cc @@ -936,6 +936,77 @@ check("NRSMat::storesubmatrix(selection) copy-on-write and size check",unchanged } +//single precision: every method with the float defaults (opt_precision) on float test problems +template +static float runf(O &o, const OptProblem *p, NRVec &x, H &h, bool *converged) +{ +h.my_verbose=0; +x.copyonwrite(); +optresult r=o.minimize(p,x,&h,NULL); +*converged=r.converged; +return r.f; +} + +static NRVec startpointf(const int n, const int seed) +{ +NRVec xd=startpoint(n,seed); +NRVec x(n); +for(int i=0; i class O, template class H> +static void compare_float(const char *name, const OptProblem *pd, const OptProblem *pf, const int n, const int seed, const int maxstep, void (*modify)(H &, H &)=NULL) +{ +NRVec xd=startpoint(n,seed); NRVec xf(n); for(int i=0; i hd; H hf; hd.my_maxstep=hf.my_maxstep=maxstep; if(modify) modify(hd,hf); +O od; O of; +bool convd, convf; +hd.my_verbose=0; xd.copyonwrite(); optresult rd=od.minimize(pd,xd,&hd,NULL); convd=rd.converged; +float ff=runf(of,pf,xf,hf,&convf); +string s=string("float ")+name+" agrees with double"; +check(s.c_str(),convd && convf && abs(ff-rd.f)<1e-5*max(1.,abs(rd.f)),ff-rd.f); +} +static void dfp(BFGS_hyperparameters &hd, BFGS_hyperparameters &hf) {hd.my_update=hf.my_update=QN_DFP; hd.my_c2=hf.my_c2=0.1;} +static void sr1(BFGS_hyperparameters &hd, BFGS_hyperparameters &hf) {hd.my_update=hf.my_update=QN_SR1;} + +static void test_float() +{ +cout <<"--- single precision (float)\n"; +testproblem tp; //f=0 at x=0, used for the derivative-free methods; its large constant term n(n-1)/2 cancels in f, so f is not + //resolved below about 3e-6 in float and gradient methods cannot reach g_thr there: they use testproblem_vqe +testproblem_vqe vd; testproblem_vqe vf; +compare_float("BFGS on testproblem_vqe",&vd,&vf,8,3,5000); +compare_float("DFP on testproblem_vqe",&vd,&vf,8,3,5000,&dfp); +compare_float("SR1 on testproblem_vqe",&vd,&vf,8,3,5000,&sr1); +compare_float("L-BFGS on testproblem_vqe",&vd,&vf,8,3,5000); +compare_float("CG on testproblem_vqe",&vd,&vf,8,3,5000); +compare_float("GDIIS on testproblem_vqe",&vd,&vf,8,3,5000); +{NelderMead_hyperparameters h; h.my_maxstep=100000; NelderMead_Optimizer o; NRVec x=startpointf(4,1); bool conv; float f=runf(o,&tp,x,h,&conv); check("float Nelder-Mead converges on testproblem n=4 (default f_thr, x_thr)",conv && f<1e-5,f);} +{Powell_hyperparameters h; Powell_Optimizer o; NRVec x=startpointf(10,1); bool conv; float f=runf(o,&tp,x,h,&conv); check("float Powell converges on testproblem n=10 (default f_thr, x_thr)",conv && f<1e-5,f);} +{Cobyla_hyperparameters h; Cobyla_Optimizer o; NRVec x=startpointf(4,1); bool conv; float f=runf(o,&tp,x,h,&conv); check("float COBYLA converges on testproblem n=4 (default rhoend)",conv && f<1e-5,f);} +{Adam_hyperparameters h; h.my_maxstep=20000; h.my_lr=1e-2; Adam_Optimizer o; NRVec x=startpointf(10,1); bool conv; float f=runf(o,&tp,x,h,&conv); check("float Adam converges on testproblem n=10",conv && f<1e-5,f);} +{//finite-difference gradient with the float default step +FiniteDifference fd(&tp); NRVec x=startpointf(10,2), g(10), gfd(10); tp.funcgrad(x,&g); fd.funcgrad(x,&gfd); gfd-=g; +check("float FiniteDifference gradient (default step 5e-3)",gfd.norm()<1e-3*g.norm(),gfd.norm()/g.norm()); +FiniteDifference fdd(&vd); FiniteDifference fdf(&vf); +compare_float("BFGS with finite-difference gradients (default steps)",&fdd,&fdf,8,3,5000); +} +{//Hessian-based methods on the VQE-like test problem, compared with double +NRVec xd=startpoint(6,3); NRVec xf(6); for(int i=0;i<6;++i) xf[i]=xd[i]; +TrustRegion_hyperparameters hd; TrustRegion_Optimizer od; double fd=run(od,&vd,xd,hd); +TrustRegion_hyperparameters hf; TrustRegion_Optimizer of; bool conv; float ff=runf(of,&vf,xf,hf,&conv); +check("float trust region on testproblem_vqe agrees with double",conv && abs(ff-fd)<1e-5*max(1.,abs(fd)),ff-fd); +} +{//least squares: NIST Misra1a in float needs Marquardt scaling (unscaled b1 steps fall below float resolution) +testproblem_misra1a p; NRVec b(2); for(int i=0;i<2;++i) b[i]=p.start(0,i); +LM_hyperparameters h; h.my_maxstep=100000; h.my_scaling=true; LM_Optimizer o; bool conv; float f=runf(o,&p,b,h,&conv); +float perr=0; for(int i=0;i<2;++i) perr=max(perr,(float)abs(b[i]-p.certified(i))/abs(p.certified(i))); +check("float Levenberg-Marquardt fits NIST Misra1a",conv && perr<1e-3,perr); +} +} + int main(int argc, char **argv) { test_problems_and_capabilities(); @@ -949,6 +1020,7 @@ test_leastsquares(); test_maxfunc_and_parameters(); test_diis(); test_copyonwrite(); +test_float(); cout <<"=== "<