diff --git a/optimize.cc b/optimize.cc index 3dcd974..cdd3cfc 100644 --- a/optimize.cc +++ b/optimize.cc @@ -17,6 +17,7 @@ */ #include "optimize.h" +#include "la_random.h" namespace LA { @@ -203,7 +204,7 @@ return true; template -bool optstate::stop(const T ff, const T gn) +bool optstate::stop(const T ff, const T gn, const int needed) { f=ff; gnorm=gn; if(hp.my_verbose) @@ -214,7 +215,7 @@ if(hp.my_verbose) } if(hp.converged(*this)) {converged=true; return true;} if(iter>=hp.my_maxstep) {if(hp.my_verbose) std::cout < +optresult firstorder_minimize(const char *name, firstorder_step &stepper, const OptProblem *o, NRVec &x, const firstorder_hyperparameters &hp, hyperparameters *achieved, gradient_estimator *est) +{ +gradient_estimator exact; +gradient_estimator &grad = est ? *est : exact; +optevaluator e(o,hp.my_maxfunc); +optstate st(name,hp,e); +const int n=x.size(); +const T lr0=grad.calibrate(e,x,hp); +if(hp.my_verbose && lr0!=hp.my_lr) std::cout < g(n); +T f=grad(e,x,g,0); + +for(;; ++st.iter) + { + if(st.stop(f,grad.exact() ? g.norm() : (T)-1,grad.cost())) break; + T lr=lr0/std::pow(1+hp.my_lr_decay*st.iter,hp.my_lr_power); + NRVec dx = stepper.step(g,st.iter+1); + dx *= -lr; + if(hp.my_weight_decay!=0) dx.axpy(-lr*hp.my_weight_decay,x); + x += dx; + T fold=f; + NRVec gnew(n); //fresh vector, the stepper may keep the previous gradient + f=grad(e,x,gnew,st.iter+1); + g=gnew; + if(grad.exact()) st.step(f-fold,dx.norm()); //for estimates, small changes are no sign of convergence + } + +if(hp.my_return_best && e.fbest stepper(hp); +return firstorder_minimize(hp.my_nesterov ? "SGD-Nesterov" : "SGD",stepper,o,x,hp,achieved); +} + + +//elementwise functions for Adam +template +static T elementmax(const T &a, const T &b) {return a>b ? a : b;} + +//m/(sqrt(v)/sqrtc2 + eps) with the bias correction sqrtc2 = sqrt(1-beta2^k) +template +class adamratio { +public: +T sqrtc2, eps; +adamratio(const T s, const T e) : sqrtc2(s), eps(e) {}; +T operator()(const T &m, const T &v) const {return m/(std::sqrt(v)/sqrtc2+eps);} +}; + + +//Adam, AdamW, AMSGrad +template +class Adam_stepper : public firstorder_step { +public: +const Adam_hyperparameters &hp; +NRVec m; //first moment estimate +NRVec v; //second moment estimate +NRVec vmax; //maximum of the second moments for AMSGrad +T beta1k, beta2k; //beta1^k, beta2^k +Adam_stepper(const Adam_hyperparameters &h) : hp(h), beta1k(1), beta2k(1) {}; +NRVec step(const NRVec &g, const int k) + { + NRVec g2(g); + g2 *= g; + if(k==1) {m = g*(1-hp.my_beta1); v = g2*(1-hp.my_beta2);} + else + { + m *= hp.my_beta1; m.axpy(1-hp.my_beta1,g); + v *= hp.my_beta2; v.axpy(1-hp.my_beta2,g2); + } + beta1k *= hp.my_beta1; + beta2k *= hp.my_beta2; + if(hp.my_amsgrad) {if(k==1) vmax=v; else vmax.call_on_me(elementmax,v);} + NRVec d = m*(1/(1-beta1k)); + d.call_on_me(adamratio(std::sqrt(1-beta2k),hp.my_eps), hp.my_amsgrad ? vmax : v); + return d; + } +}; + + +template +optresult Adam_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +Adam_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("Adam cannot search for saddle points"); +Adam_stepper stepper(hp); +return firstorder_minimize(hp.my_amsgrad ? "AMSGrad" : (hp.my_weight_decay!=0 ? "AdamW" : "Adam"),stepper,o,x,hp,achieved); +} + + + +//SPSA gradient estimate (g_hat)_i = (f(x+c_k Delta) - f(x-c_k Delta))/(2 c_k Delta_i), Delta_i = +-1, averaged over resamplings +template +class SPSA_estimator : public gradient_estimator { +public: +const SPSA_hyperparameters &hp; +SPSA_estimator(const SPSA_hyperparameters &h) : hp(h) {}; +T operator()(optevaluator &e, const NRVec &x, NRVec &g, const int k); +int cost() const {return 2*hp.my_resamplings;} +bool exact() const {return false;} +T calibrate(optevaluator &e, const NRVec &x, const firstorder_hyperparameters &h); +}; + +//returns the average of f(x+c_k Delta) and f(x-c_k Delta) as an estimate of f(x) +template +T SPSA_estimator::operator()(optevaluator &e, const NRVec &x, NRVec &g, const int k) +{ +const int n=x.size(); +const T ck=hp.my_c/std::pow((T)(k+1),hp.my_gamma); +g.clear(); +T fsum=0; +for(int r=0; r delta(n); //fresh vector, element access is safe + for(int i=0; i xp(x), xm(x); + xp.axpy(ck,delta); + xm.axpy(-ck,delta); + T fp=e(xp); + T fm=e(xm); + g.axpy((fp-fm)/(2*ck),delta); //1/Delta_i = Delta_i + fsum+=fp+fm; + } +g /= (T)hp.my_resamplings; +return fsum/(2*hp.my_resamplings); +} + +//learning rate such that the first step lr*|g_hat| has on average the target length +template +T SPSA_estimator::calibrate(optevaluator &e, const NRVec &x, const firstorder_hyperparameters &h) +{ +if(hp.my_calibrate<=0) return hp.my_lr; +T gsum=0; +int ncal=0; +for(; ncal=2*cost(); ++ncal) //leave the budget for the first estimate + { + NRVec g(x.size()); + (*this)(e,x,g,0); + gsum+=g.norm(); + } +if(ncal==0 || gsum==0) return hp.my_lr; +return hp.my_target_step*ncal/gsum; +} + + +//plain gradient step for SPSA +template +class plain_stepper : public firstorder_step { +public: +NRVec step(const NRVec &g, const int k) {return g;} +}; + + +template +optresult SPSA_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +SPSA_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("SPSA cannot search for saddle points"); +if(hp.my_resamplings<1) laerror("SPSA needs resamplings >= 1"); +plain_stepper stepper; +SPSA_estimator est(hp); +return firstorder_minimize("SPSA",stepper,o,x,hp,achieved,&est); +} + + + //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; @@ -750,6 +954,9 @@ template class LBFGS_Optimizer; template class CG_Optimizer; template class NelderMead_Optimizer; template class Powell_Optimizer; +template class SGD_Optimizer; +template class Adam_Optimizer; +template class SPSA_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 6c854dc..cfd6d3f 100644 --- a/optimize.h +++ b/optimize.h @@ -145,7 +145,7 @@ T dx; //|x change| in the last step bool stepdone; //df and dx are valid bool converged; optstate(const char *nm, const hyperparameters &h, const optevaluator &ee) : name(nm), hp(h), e(ee), iter(0), f(0), gnorm(-1), df(0), dx(0), stepdone(false), converged(false) {}; -bool stop(const T ff, const T gn= -1); //to be called at the beginning of each iteration, returns true if the iterations should finish +bool stop(const T ff, const T gn= -1, const int needed=1); //to be called at the beginning of each iteration, returns true if the iterations should finish; needed = no. of function evaluations of the next iteration void step(const T dff, const T dxx) {df=std::abs(dff); dx=dxx; stepdone=true;} //record a step done optresult finish(hyperparameters *achieved) const; //final report, fills achieved if not NULL }; @@ -197,6 +197,44 @@ template optresult linesearch_minimize(const char *name, linesearch_direction &dir, const OptProblem *o, NRVec &x, const linesearch_hyperparameters &hp, hyperparameters *achieved); +//parameters of first-order methods with a learning rate instead of line search, suitable also for noisy (stochastic) functions and gradients +template +class firstorder_hyperparameters : public gradient_hyperparameters { +public: +T my_lr; //learning rate +T my_lr_decay; //learning rate in iteration k is lr/(1+lr_decay*k)^lr_power +T my_lr_power; +T my_weight_decay; //decoupled weight decay x -= lr*weight_decay*x (Loshchilov & Hutter, ICLR 2019) +bool my_return_best; //return the point with the lowest function value found instead of the last one; only for (nearly) deterministic f, the lowest of noisy values is biased by the noise +firstorder_hyperparameters() : my_lr(1e-3), my_lr_decay(0), my_lr_power(1), my_weight_decay(0), my_return_best(false) {this->my_maxstep=1000;}; +virtual void setfrom(const hyperparameters &p) {const firstorder_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else gradient_hyperparameters::setfrom(p);} +}; + +//update rule of a first-order method, used by the common driver firstorder_minimize() +template +class firstorder_step { +public: +virtual NRVec step(const NRVec &g, const int k) =0; //direction d of the step x -= lr*d in iteration k=1,2,... with gradient g +virtual ~firstorder_step() {}; +}; + +//function value and gradient for firstorder_minimize(), the default is the exact evaluation, derived classes estimate them (SPSA) +template +class gradient_estimator { +public: +virtual T operator()(optevaluator &e, const NRVec &x, NRVec &g, const int k) {return e(x,&g);} //returns f at x and the gradient g, k=0,1,... is the iteration +virtual int cost() const {return 1;} //no. of function evaluations of one call +virtual bool exact() const {return true;} //with estimates no convergence test is done, small gradients or changes can occur by chance +virtual T calibrate(optevaluator &e, const NRVec &x, const firstorder_hyperparameters &hp) {return hp.my_lr;} //learning rate to be used +virtual ~gradient_estimator() {}; +}; + +//common driver of first-order methods; est==NULL means exact gradients +//the learning rate actually used is returned in achieved, if it is of firstorder_hyperparameters type +template +optresult firstorder_minimize(const char *name, firstorder_step &stepper, const OptProblem *o, NRVec &x, const firstorder_hyperparameters &hp, hyperparameters *achieved, gradient_estimator *est=NULL); + + //specific methods //parameters common to quasi-Newton methods @@ -311,6 +349,69 @@ virtual ~Powell_Optimizer() {}; }; +//stochastic gradient descent with momentum v = momentum*v + g, step d = v or d = g + momentum*v (Nesterov) (Sutskever et al., ICML 2013) +template +class SGD_hyperparameters : public firstorder_hyperparameters { +public: +T my_momentum; +bool my_nesterov; +SGD_hyperparameters() : my_momentum(0.9), my_nesterov(false) {this->my_lr=1e-2;}; +virtual void setfrom(const hyperparameters &p) {const SGD_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else firstorder_hyperparameters::setfrom(p);} +}; + +template +class SGD_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~SGD_Optimizer() {}; +}; + + +//Adam (Kingma & Ba, ICLR 2015), with weight_decay>0 AdamW, optionally AMSGrad (Reddi et al., ICLR 2018); beta1=0 gives RMSprop +template +class Adam_hyperparameters : public firstorder_hyperparameters { +public: +T my_beta1; //decay rate of the first moment estimate +T my_beta2; //decay rate of the second moment estimate +T my_eps; //regularization of the division by the second moment +bool my_amsgrad; //use the maximum of the past second moments +Adam_hyperparameters() : my_beta1(0.9), my_beta2(0.999), my_eps(1e-8), my_amsgrad(false) {}; +virtual void setfrom(const hyperparameters &p) {const Adam_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else firstorder_hyperparameters::setfrom(p);} +}; + +template +class Adam_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~Adam_Optimizer() {}; +}; + + +//simultaneous perturbation stochastic approximation (Spall, IEEE Trans. Autom. Control 37, 332 (1992); IEEE Trans. Aerosp. Electron. Syst. 34, 817 (1998)) +//derivative-free, 2 function evaluations per gradient estimate independently of the dimension, robust to noise (VQE) +//gain sequences: a_k = lr/(1+lr_decay*k)^lr_power, i.e. Spall's a/(k+1+A)^alpha with lr=a/(1+A)^alpha, lr_decay=1/(1+A), lr_power=alpha; c_k = c/(k+1)^gamma +//the gradient and f are estimates, therefore no convergence test is done and the iterations run until maxstep or maxfunc +//the first estimate at the starting point is always done, even if maxfunc is smaller than its 2*resamplings evaluations +template +class SPSA_hyperparameters : public firstorder_hyperparameters { +public: +T my_c; //perturbation size, of the order of the noise of f +T my_gamma; //decay exponent of the perturbation size +int my_resamplings; //no. of gradient estimates averaged in one iteration +int my_calibrate; //no. of gradient estimates at the starting point for the calibration of lr, 0 = no calibration +T my_target_step; //calibration: lr is set so that the first step has this length on average +SPSA_hyperparameters() : my_c(0.1), my_gamma(0.101), my_resamplings(1), my_calibrate(20), my_target_step(0.1) {this->my_lr_decay=0.01; this->my_lr_power=0.602; this->my_g_thr=0;}; +virtual void setfrom(const hyperparameters &p) {const SPSA_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else firstorder_hyperparameters::setfrom(p);} +}; + +template +class SPSA_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~SPSA_Optimizer() {}; +}; + + }//namespace #endif /* _OPTIMIZE_H_ */ diff --git a/t.cc b/t.cc index b3b5f83..e5585f0 100644 --- a/t.cc +++ b/t.cc @@ -4847,11 +4847,16 @@ testproblem problem; //LBFGS_hyperparameters param; //tested //CG_hyperparameters param; //tested //NelderMead_hyperparameters param; //tested -Powell_hyperparameters param; +//Powell_hyperparameters param; //tested +//SGD_hyperparameters param; //tested +//Adam_hyperparameters param; //tested +SPSA_hyperparameters param; param.my_maxstep=100; param.my_f_thr=1e-6; param.my_x_thr=1e-3; -//param.my_g_thr=1e-3; //gradient methods +param.my_g_thr=1e-3; //gradient methods +param.my_lr=1e-2; //first-order methods +param.my_maxstep=2000; //first-order methods //param.my_update=QN_BFGS; //tested //param.my_update=QN_DFP; //tested //param.my_update=QN_SR1; //tested @@ -4860,13 +4865,16 @@ cout << "initial " < o; //tested //CG_Optimizer o; //tested //NelderMead_Optimizer o; //tested -Powell_Optimizer o; -//gradient_hyperparameters achieved; //gradient methods -hyperparameters achieved; +//Powell_Optimizer o; //tested +//SGD_Optimizer o; //tested +//Adam_Optimizer o; //tested +SPSA_Optimizer o; +gradient_hyperparameters achieved; //gradient methods +//hyperparameters achieved; //derivative-free methods optresult res=o.minimize(&problem,x,¶m,&achieved); cout << "optimized " < + NRVec& call_on_me(F _F){ + NOT_GPU(*this); + + copyonwrite(); + for(int i=0; i + NRVec& call_on_me(F _F, const NRVec &rhs){ +#ifdef DEBUG + if(nn!=rhs.nn) laerror("incompatible dimensions in call_on_me"); +#endif + NOT_GPU(*this); + NOT_GPU(rhs); + + copyonwrite(); + for(int i=0; i