From 3b7be8bc43b0070e257a74ee6a7e237f2c915ba0 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Mon, 26 Oct 2026 12:39:09 +0100 Subject: [PATCH] optimize: more gradient-based methods vibe-coded --- optimize.cc | 262 ++++++++++++++++++++++++++++++++++++++++++++-------- optimize.h | 79 +++++++++++++++- t.cc | 13 ++- 3 files changed, 308 insertions(+), 46 deletions(-) diff --git a/optimize.cc b/optimize.cc index c1b8053..1532c2a 100644 --- a/optimize.cc +++ b/optimize.cc @@ -151,18 +151,54 @@ return 0; } -//BFGS update of the inverse Hessian approximation H, cf. Nocedal & Wright eq. 6.17 -//H+ = (I - rho s y^T) H (I - rho y s^T) + rho s s^T, rho=1/(y^T s) -//expanded: H+ = H + (rho^2 y^T H y + rho) s s^T - rho (s (Hy)^T + Hy s^T) = H + s u^T + u s^T with u = (rho^2 y^T H y + rho)/2 s - rho Hy +//quasi-Newton update of the inverse Hessian approximation H, returns false if the update was skipped +//Hinit = H is still the initial guess, it is then scaled by s.y/y.y (Nocedal & Wright eq. 6.20) +//for SR1 the scaling makes (s-Hy).y=0, so its first update is skipped and the next one starts from the scaled H +//rank-1 terms are done via rank2update(a,x,x) = 2a x x^T template -static void BFGS_update(NRSMat &H, const NRVec &s, const NRVec &y, const T sy) +static bool QN_update_H(NRSMat &H, bool &Hinit, const NRVec &s, const NRVec &y, const BFGS_hyperparameters &hp) { -T rho = 1/sy; -NRVec Hy = H*y; -T yHy = y*Hy; -NRVec u = s*((rho*rho*yHy + rho)/2); -u.axpy(-rho,Hy); -H.rank2update((T)1,s,u); +const T sy = s*y; +const T snorm=s.norm(); +const T ynorm=y.norm(); +const bool curvature = sy > hp.my_sy_thr*snorm*ynorm; +if(Hinit && curvature) {H = (T)(sy/(ynorm*ynorm)); Hinit=false;} +switch(hp.my_update) + { + case QN_BFGS: //H+ = (I - rho s y^T) H (I - rho y s^T) + rho s s^T = H + s u^T + u s^T, rho=1/(y^T s), u = (rho^2 y^T H y + rho)/2 s - rho Hy + { + if(!curvature) return false; + T rho = 1/sy; + NRVec Hy = H*y; + NRVec u = s*((rho*rho*(y*Hy) + rho)/2); + u.axpy(-rho,Hy); + H.rank2update((T)1,s,u); + } + break; + case QN_DFP: //H+ = H - Hy (Hy)^T/(y^T H y) + s s^T/(y^T s) + { + if(!curvature) return false; + NRVec Hy = H*y; + T yHy = y*Hy; + if(yHy<=0) return false; + H.rank2update(1/(2*sy),s,s); + H.rank2update(-1/(2*yHy),Hy,Hy); + } + break; + case QN_SR1: //H+ = H + w w^T/(w^T y), w = s - Hy; H may become indefinite + { + NRVec w = H*y; + w *= (T)-1; + w += s; + T wy = w*y; + if(std::abs(wy) <= hp.my_sr1_thr*w.norm()*ynorm) return false; + H.rank2update(1/(2*wy),w,w); + } + break; + default: laerror("unknown quasi-Newton update"); + } +Hinit=false; +return true; } @@ -198,73 +234,221 @@ return optresult(f,converged); template -optresult BFGS_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +optresult linesearch_minimize(const char *name, linesearch_direction &dir, const OptProblem *o, NRVec &x, const linesearch_hyperparameters &hp, hyperparameters *achieved) { -BFGS_hyperparameters hp; -if(par) hp.setfrom(*par); -if(hp.my_saddle_order) laerror("BFGS cannot search for saddle points"); const int verbose=hp.my_verbose; optevaluator e(o,hp.my_maxfunc); -optstate st("BFGS",hp,e); +optstate st(name,hp,e); const int n=x.size(); NRVec g(n); T f=e(x,&g); -NRSMat H(n); -H=1.; -bool Hinit=true; //H is the initial guess, rescale it after the first step +T alpha=0, dphi=0; //of the previous step for(;; ++st.iter) { if(st.stop(f,g.norm())) break; - NRVec p = H*g; - p *= (T)-1; - if(p*g >= 0) //H lost positive definiteness numerically, restart + NRVec p = dir.direction(g); + bool sd=dir.steepest(); //steepest descent step + if(p*g >= 0) { - if(verbose) std::cout <<"BFGS: not a descent direction, resetting inverse Hessian\n"; - H=1.; Hinit=true; + if(dir.indefinite()) {if(verbose) std::cout <0) alpha1=alpha*dphi/dphi_new; + else alpha1= sd ? 1/p.norm() : (T)1; + dphi=dphi_new; + NRVec xold(x), gold(g); T fold=f; - //on the first step (or after restart) try unit step length in x, otherwise the Newton step - T alpha=linesearch_wolfe(e,x,f,g,p,Hinit ? 1/p.norm() : (T)1,hp); + alpha=linesearch_wolfe(e,x,f,g,p,alpha1,hp); if(alpha==0) { if(e.exhausted()) continue; //reported by stop() - if(Hinit) {if(verbose) std::cout <<"BFGS: line search failed along steepest descent, giving up\n"; break;} - if(verbose) std::cout <<"BFGS: line search failed, resetting inverse Hessian\n"; - H=1.; Hinit=true; + if(sd) {if(verbose) std::cout < s = x-xold; NRVec y = g-gold; - T sy = s*y; - T snorm=s.norm(); - T ynorm=y.norm(); - if(verbose>1) std::cout <<"BFGS step alpha= "<1) std::cout < +class QN_direction : public linesearch_direction { +public: +const BFGS_hyperparameters &hp; +NRSMat H; +bool Hinit; //H is the initial guess, rescale it after the first step +QN_direction(const BFGS_hyperparameters &h, const int n) : hp(h), H(n), Hinit(true) {H=1.;}; +NRVec direction(const NRVec &g) {NRVec p = H*g; p *= (T)-1; return p;} +bool update(const NRVec &s, const NRVec &y, const NRVec &p, const NRVec &g) {return QN_update_H(H,Hinit,s,y,hp);} +void reset() {H=1.; Hinit=true;} +bool steepest() const {return Hinit;} +bool indefinite() const {return hp.my_update==QN_SR1;} +}; + + +template +optresult BFGS_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +static const char *names[]={"BFGS","DFP","SR1"}; +BFGS_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("quasi-Newton with line search cannot search for saddle points"); +QN_direction dir(hp,x.size()); +return linesearch_minimize(names[hp.my_update],dir,o,x,hp,achieved); +} + + +//L-BFGS: last m pairs s = x_{k+1}-x_k, y = g_{k+1}-g_k in a ring buffer +template +class LBFGS_direction : public linesearch_direction { +public: +const LBFGS_hyperparameters &hp; +int m; //capacity +int count; //no. of pairs stored +int newest; //index of the newest pair +NRVec > S,Y; +NRVec rho; //1/(s.y) +T gamma; //initial inverse Hessian H0 = gamma I, gamma = s.y/y.y of the newest pair (Nocedal & Wright eq. 7.20) +LBFGS_direction(const LBFGS_hyperparameters &h) : hp(h), m(h.my_memory), count(0), newest(-1), S(m), Y(m), rho(m), gamma(1) {}; +void push(const NRVec &s, const NRVec &y, const T sy) {newest=(newest+1)%m; S[newest]=s; Y[newest]=y; rho[newest]=1/sy; gamma=sy/(y*y); if(count direction(const NRVec &g); +bool update(const NRVec &s, const NRVec &y, const NRVec &p, const NRVec &g) + { + T sy=s*y; + if(sy <= hp.my_sy_thr*s.norm()*y.norm()) return false; //curvature condition violated numerically + push(s,y,sy); + return true; + } +void reset() {count=0; newest= -1; gamma=1;} +bool steepest() const {return count==0;} +}; + +//returns -H g by the two-loop recursion (Nocedal & Wright alg. 7.4) +template +NRVec LBFGS_direction::direction(const NRVec &g) +{ +NRVec q(g); +NRVec alpha(count); +int k=newest; +for(int i=0; i=0; --i) //oldest to newest + { + T beta = rho[k]*(Y[k]*q); + q.axpy(alpha[i]-beta,S[k]); + k=(k+1)%m; + } +q *= (T)-1; +return q; +} + + +template +optresult LBFGS_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +LBFGS_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("L-BFGS cannot search for saddle points"); +if(hp.my_memory<1) laerror("L-BFGS needs memory >= 1"); +LBFGS_direction dir(hp); +return linesearch_minimize("L-BFGS",dir,o,x,hp,achieved); +} + + +//nonlinear CG: p_{k+1} = -g_{k+1} + beta p_k +template +class CG_direction : public linesearch_direction { +public: +const CG_hyperparameters &hp; +int restart_interval; +int count; //no. of steps since the last restart +bool restart; //next direction is steepest descent +T beta; +T gg; //g.g of the gradient of the last direction +NRVec d; //last search direction +CG_direction(const CG_hyperparameters &h, const int n) : hp(h), restart_interval(h.my_restart>0 ? h.my_restart : n), count(0), restart(true), beta(0), gg(0) {}; +NRVec direction(const NRVec &g) {NRVec p = g*(T)-1; if(!restart) p.axpy(beta,d); gg=g*g; return p;} +bool update(const NRVec &s, const NRVec &y, const NRVec &p, const NRVec &g); +void reset() {restart=true; count=0;} +bool steepest() const {return restart;} +bool alpha_from_previous() const {return true;} +}; + +template +bool CG_direction::update(const NRVec &s, const NRVec &y, const NRVec &p, const NRVec &g) +{ +d=p; +T gg1=g*g; +T gy=g*y; +T dy=p*y; +restart= ++count>=restart_interval || (hp.my_powell_thr>0 && std::abs(gg1-gy) >= hp.my_powell_thr*gg1); //g_{k+1}.g_k = gg1-gy +bool curvature= dy>0; +if(!restart) switch(hp.my_beta) + { + case CG_FR: beta=gg1/gg; break; + case CG_PR: beta=std::max((T)0,gy/gg); break; + case CG_HS: if(curvature) beta=std::max((T)0,gy/dy); else restart=true; break; + case CG_DY: if(curvature) beta=gg1/dy; else restart=true; break; + case CG_HZ: if(curvature) beta=std::max((gy-2*(p*g)*(y*y)/dy)/dy, -1/(p.norm()*std::min((T)0.01,std::sqrt(gg)))); else restart=true; break; + default: laerror("unknown CG beta formula"); + } +if(restart) + { + count=0; + if(hp.my_verbose>1) std::cout <<"CG restart\n"; + } +return true; +} + + +template +optresult CG_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +static const char *names[]={"CG-FR","CG-PR+","CG-HS+","CG-DY","CG-HZ"}; +CG_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("CG cannot search for saddle points"); +CG_direction dir(hp,x.size()); +return linesearch_minimize(names[hp.my_beta],dir,o,x,hp,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; template class testproblem; template class BFGS_Optimizer; +template class LBFGS_Optimizer; +template class CG_Optimizer; }//namespace diff --git a/optimize.h b/optimize.h index da3fbef..76d57fe 100644 --- a/optimize.h +++ b/optimize.h @@ -172,14 +172,46 @@ template T linesearch_wolfe(optevaluator &e, NRVec &x, T &f, NRVec &g, const NRVec &p, const T alpha1, const linesearch_hyperparameters &lp); +//search direction of a line-search method, used by the common driver linesearch_minimize() +template +class linesearch_direction { +public: +virtual NRVec direction(const NRVec &g) =0; //search direction at the point with gradient g +virtual bool update(const NRVec &s, const NRVec &y, const NRVec &p, const NRVec &g) =0; //after the step s=alpha*p with gradient change y and new gradient g; returns false if the update was skipped +virtual void reset() =0; //forget the accumulated information +virtual bool steepest() const =0; //no information accumulated, direction() gives steepest descent +virtual bool indefinite() const {return false;} //a non-descent direction is legitimate (SR1), take steepest descent step without reset +virtual bool alpha_from_previous() const {return false;} //initial step length from the previous one (Nocedal & Wright eq. 3.60), otherwise 1 (unit step in x for steepest descent) +virtual ~linesearch_direction() {}; +}; + +//common driver of line-search methods +template +optresult linesearch_minimize(const char *name, linesearch_direction &dir, const OptProblem *o, NRVec &x, const linesearch_hyperparameters &hp, hyperparameters *achieved); + + //specific methods +//parameters common to quasi-Newton methods template -class BFGS_hyperparameters : public linesearch_hyperparameters { +class quasinewton_hyperparameters : public linesearch_hyperparameters { public: -T my_sy_thr; //skip the update if s.y <= sy_thr*|s|*|y| (curvature condition violated numerically) -BFGS_hyperparameters() : my_sy_thr(1e-12) {}; -virtual void setfrom(const hyperparameters &p) {const BFGS_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else linesearch_hyperparameters::setfrom(p);} +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) {}; +virtual void setfrom(const hyperparameters &p) {const quasinewton_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else linesearch_hyperparameters::setfrom(p);} +}; + + +//quasi-Newton methods with line search, differing by the update of the inverse Hessian approximation +enum QN_update {QN_BFGS, QN_DFP, QN_SR1}; + +template +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) {}; +virtual void setfrom(const hyperparameters &p) {const BFGS_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else quasinewton_hyperparameters::setfrom(p);} }; template @@ -190,6 +222,45 @@ virtual ~BFGS_Optimizer() {}; }; +//limited-memory BFGS (Nocedal & Wright chap. 7.2), for large problems +template +class LBFGS_hyperparameters : public quasinewton_hyperparameters { +public: +int my_memory; //no. of stored pairs of steps and gradient changes +LBFGS_hyperparameters() : my_memory(10) {}; +virtual void setfrom(const hyperparameters &p) {const LBFGS_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else quasinewton_hyperparameters::setfrom(p);} +}; + +template +class LBFGS_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~LBFGS_Optimizer() {}; +}; + + +//nonlinear conjugate gradient methods (Nocedal & Wright chap. 5.2; Hager & Zhang, SIAM J. Optim. 16, 170 (2005)) +//differing by the formula for beta in p_{k+1} = -g_{k+1} + beta p_k; PR and HS are used in the nonnegative variants PR+ and HS+ +enum CG_beta {CG_FR, CG_PR, CG_HS, CG_DY, CG_HZ}; + +template +class CG_hyperparameters : public linesearch_hyperparameters { +public: +CG_beta my_beta; //Fletcher-Reeves, Polak-Ribiere+, Hestenes-Stiefel+, Dai-Yuan, Hager-Zhang +int my_restart; //restart with steepest descent every my_restart iterations, <=0 means the problem dimension +T my_powell_thr; //restart if |g_{k+1}.g_k| >= powell_thr*|g_{k+1}|^2 (Nocedal & Wright eq. 5.52), <=0 disables +CG_hyperparameters() : my_beta(CG_PR), my_restart(0), my_powell_thr(0.2) {this->my_c2=0.1;}; //CG needs a more accurate line search +virtual void setfrom(const hyperparameters &p) {const CG_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else linesearch_hyperparameters::setfrom(p);} +}; + +template +class CG_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~CG_Optimizer() {}; +}; + + }//namespace #endif /* _OPTIMIZE_H_ */ diff --git a/t.cc b/t.cc index 252877f..ce55aa1 100644 --- a/t.cc +++ b/t.cc @@ -4843,14 +4843,21 @@ cin >>n; NRVec x(n); x.randomize(1.); testproblem problem; -BFGS_hyperparameters param; +//BFGS_hyperparameters param; //tested +//LBFGS_hyperparameters param; //tested +CG_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_update=QN_BFGS; //tested +//param.my_update=QN_DFP; //tested +//param.my_update=QN_SR1; //tested cout << "initial " < o; -BFGS_hyperparameters achieved; +//BFGS_Optimizer o; //tested +//LBFGS_Optimizer o; //tested +CG_Optimizer o; +gradient_hyperparameters achieved; optresult res=o.minimize(&problem,x,¶m,&achieved); cout << "optimized " <