From ae7a40bd771f8eb26f83dfb77c103e769c5af8d4 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Mon, 26 Oct 2026 11:34:03 +0100
Subject: [PATCH] optimize: reorganization of hyperparameters and optevaluator
 and optstate

---
 optimize.cc | 126 ++++++++++++++++++++++++++--------------------------
 optimize.h  |  82 +++++++++++++++++++++++++++++-----
 2 files changed, 134 insertions(+), 74 deletions(-)

diff --git a/optimize.cc b/optimize.cc
index 641d3cf..c1b8053 100644
--- a/optimize.cc
+++ b/optimize.cc
@@ -70,14 +70,19 @@ return r;
 
 
 template<typename T>
-T linesearch_wolfe(const OptProblem<T> *o, NRVec<T> &x, T &f, NRVec<T> &g, const NRVec<T> &p, T alpha1, const T c1, const T c2, const T alphamax, const int maxiter, int &nfunc, const int verbose)
+T linesearch_wolfe(optevaluator<T> &e, NRVec<T> &x, T &f, NRVec<T> &g, const NRVec<T> &p, const T alpha1, const linesearch_hyperparameters<T> &lp)
 {
+const T c1=lp.my_c1;
+const T c2=lp.my_c2;
+const int verbose=lp.my_verbose;
+const int maxiter= std::min(lp.my_maxlinesearch,e.remaining());
 const int n=x.size();
 const NRVec<T> x0(x);
 const T f0=f;
 const T d0= g*p;
 if(d0>=0) {if(verbose) std::cout <<"linesearch: not a descent direction\n"; return 0;}
-if(alpha1>alphamax) alpha1=alphamax;
+const T alphamax=lp.my_maxsteplength/p.norm();
+T alpha=std::min(alpha1,alphamax);
 
 //data of trial points; vectors are always freshly allocated, so that assignments (which share data) are safe
 T alpha0=0, f_alpha0=f0, d_alpha0=d0; NRVec<T> x_alpha0(x0), g_alpha0(g); //previous trial point
@@ -88,28 +93,28 @@ int it;
 //bracketing phase
 for(it=0; it<maxiter; ++it)
 	{
-	NRVec<T> xt = x0 + p*alpha1;
+	NRVec<T> xt = x0 + p*alpha;
 	NRVec<T> gt(n);
-	T ft=o->funcgrad(xt,&gt); ++nfunc;
+	T ft=e(xt,&gt);
 	T dt= gt*p;
-	if(verbose>1) std::cout <<"linesearch bracket "<<it<<" alpha= "<<alpha1<<" f= "<<ft<<" dphi= "<<dt<<std::endl;
-	if(ft > f0 + c1*alpha1*d0 || (it>0 && ft>=f_alpha0))
+	if(verbose>1) std::cout <<"linesearch bracket "<<it<<" alpha= "<<alpha<<" f= "<<ft<<" dphi= "<<dt<<std::endl;
+	if(ft > f0 + c1*alpha*d0 || (it>0 && ft>=f_alpha0))
 		{
 		alo=alpha0; flo=f_alpha0; dlo=d_alpha0; xlo=x_alpha0; glo=g_alpha0;
-		ahi=alpha1; fhi=ft; dhi=dt;
+		ahi=alpha; fhi=ft; dhi=dt;
 		bracketed=true;
 		break;
 		}
-	if(std::abs(dt) <= -c2*d0 || alpha1>=alphamax) {x=xt; f=ft; g=gt; return alpha1;} //strong Wolfe satisfied or cannot go further
+	if(std::abs(dt) <= -c2*d0 || alpha>=alphamax) {x=xt; f=ft; g=gt; return alpha;} //strong Wolfe satisfied or cannot go further
 	if(dt>=0)
 		{
-		alo=alpha1; flo=ft; dlo=dt; xlo=xt; glo=gt;
+		alo=alpha; flo=ft; dlo=dt; xlo=xt; glo=gt;
 		ahi=alpha0; fhi=f_alpha0; dhi=d_alpha0;
 		bracketed=true;
 		break;
 		}
-	alpha0=alpha1; f_alpha0=ft; d_alpha0=dt; x_alpha0=xt; g_alpha0=gt;
-	alpha1 *= 2; if(alpha1>alphamax) alpha1=alphamax;
+	alpha0=alpha; f_alpha0=ft; d_alpha0=dt; x_alpha0=xt; g_alpha0=gt;
+	alpha *= 2; if(alpha>alphamax) alpha=alphamax;
 	}
 
 //zoom phase, invariant: alo satisfies sufficient decrease and has the lowest f found so far, dlo*(ahi-alo)<0
@@ -118,7 +123,7 @@ if(bracketed) for(++it; it<maxiter; ++it)
 	T a=cubic_interpolate(alo,flo,dlo,ahi,fhi,dhi);
 	NRVec<T> xt = x0 + p*a;
 	NRVec<T> gt(n);
-	T ft=o->funcgrad(xt,&gt); ++nfunc;
+	T ft=e(xt,&gt);
 	T dt= gt*p;
 	if(verbose>1) std::cout <<"linesearch zoom "<<it<<" alpha= "<<a<<" f= "<<ft<<" dphi= "<<dt<<std::endl;
 	if(ft > f0 + c1*a*d0 || ft>=flo)
@@ -161,38 +166,57 @@ H.rank2update((T)1,s,u);
 }
 
 
+template<typename T>
+bool optstate<T>::stop(const T ff, const T gn)
+{
+f=ff; gnorm=gn;
+if(hp.my_verbose)
+	{
+	std::cout <<name<<" iter "<<iter<<" f= "<<f;
+	if(gnorm>=0) std::cout <<" |g|= "<<gnorm;
+	std::cout <<std::endl;
+	}
+if(hp.converged(*this)) {converged=true; return true;}
+if(iter>=hp.my_maxstep) {if(hp.my_verbose) std::cout <<name<<": maximum number of iterations reached\n"; return true;}
+if(e.exhausted()) {if(hp.my_verbose) std::cout <<name<<": maximum number of function evaluations reached\n"; return true;}
+return false;
+}
+
+
+template<typename T>
+optresult<T> optstate<T>::finish(hyperparameters<T> *achieved) const
+{
+if(hp.my_verbose)
+	{
+	std::cout <<name<<(converged?" converged in ":" NOT converged after ")<<iter<<" iterations, "<<e.nfunc<<" function evaluations, f= "<<f;
+	if(gnorm>=0) std::cout <<" |g|= "<<gnorm;
+	std::cout <<std::endl;
+	}
+if(achieved) {achieved->setfrom(hp); achieved->setreached(*this);}
+return optresult<T>(f,converged);
+}
+
+
 template<typename T>
 optresult<T> BFGS_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
 {
-//collect parameters, those not provided by the caller are taken from defaults
 BFGS_hyperparameters<T> hp;
-if(par)
-	{
-	static_cast<hyperparameters<T> &>(hp) = *par;
-	const gradient_hyperparameters<T> *gpar = dynamic_cast<const gradient_hyperparameters<T> *>(par);
-	if(gpar) static_cast<gradient_hyperparameters<T> &>(hp) = *gpar;
-	const BFGS_hyperparameters<T> *bpar = dynamic_cast<const BFGS_hyperparameters<T> *>(par);
-	if(bpar) hp = *bpar;
-	}
+if(par) hp.setfrom(*par);
+if(hp.my_saddle_order) laerror("BFGS cannot search for saddle points");
 const int verbose=hp.my_verbose;
+optevaluator<T> e(o,hp.my_maxfunc);
+optstate<T> st("BFGS",hp,e);
 
 const int n=x.size();
 NRVec<T> g(n);
-T f=o->funcgrad(x,&g);
-int nfunc=1;
+T f=e(x,&g);
 NRSMat<T> H(n);
 H=1.;
 bool Hinit=true; //H is the initial guess, rescale it after the first step
-bool converged=false;
-T last_df=0, last_dx=0;
 
-int iter;
-for(iter=0; iter<hp.my_maxstep; ++iter)
+for(;; ++st.iter)
 	{
-	T gnorm=g.norm();
-	if(verbose) std::cout <<"BFGS iter "<<iter<<" f= "<<f<<" |g|= "<<gnorm<<std::endl;
-	if(gnorm<hp.my_g_thr) {converged=true; break;}
-	if(hp.my_maxfunc>0 && nfunc>=hp.my_maxfunc) {if(verbose) std::cout <<"BFGS: maximum number of function evaluations reached\n"; break;}
+	if(st.stop(f,g.norm())) break;
 
 	NRVec<T> p = H*g;
 	p *= (T)-1;
@@ -205,15 +229,11 @@ for(iter=0; iter<hp.my_maxstep; ++iter)
 
 	NRVec<T> xold(x), gold(g);
 	T fold=f;
-	T pnorm=p.norm();
-	T alpha1 = Hinit ? 1/pnorm : 1; //on the first step (or after restart) try unit step length in x, otherwise the Newton step
-	T alphamax = hp.my_maxsteplength/pnorm;
-	int maxls=hp.my_maxlinesearch;
-	if(hp.my_maxfunc>0 && hp.my_maxfunc-nfunc<maxls) maxls=hp.my_maxfunc-nfunc;
-	T alpha=linesearch_wolfe(o,x,f,g,p,alpha1,hp.my_c1,hp.my_c2,alphamax,maxls,nfunc,verbose);
+	//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);
 	if(alpha==0)
 		{
-		if(hp.my_maxfunc>0 && nfunc>=hp.my_maxfunc) continue; //reported at the top of the loop
+		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;
@@ -232,39 +252,17 @@ for(iter=0; iter<hp.my_maxstep; ++iter)
 		BFGS_update(H,s,y,sy);
 		}
 	else if(verbose) std::cout <<"BFGS: curvature condition violated, skipping update\n";
-
-	last_df=std::abs(f-fold); last_dx=snorm;
-	if(last_df<hp.my_f_thr && last_dx<hp.my_x_thr) {converged=true; ++iter; break;}
+	st.step(f-fold,snorm);
 	}
 
-if(verbose)
-	{
-	if(converged) std::cout <<"BFGS converged in "<<iter<<" iterations, "<<nfunc<<" function evaluations, f= "<<f<<" |g|= "<<g.norm()<<std::endl;
-	else std::cout <<"BFGS NOT converged after "<<iter<<" iterations, "<<nfunc<<" function evaluations, f= "<<f<<" |g|= "<<g.norm()<<std::endl;
-	}
-
-if(achieved)
-	{
-	BFGS_hyperparameters<T> ach(hp);
-	ach.my_maxstep=iter;
-	ach.my_maxfunc=nfunc;
-	ach.my_f_thr=last_df;
-	ach.my_x_thr=last_dx;
-	ach.my_g_thr=g.norm();
-	*achieved = static_cast<hyperparameters<T> &>(ach);
-	gradient_hyperparameters<T> *gach = dynamic_cast<gradient_hyperparameters<T> *>(achieved);
-	if(gach) *gach = static_cast<gradient_hyperparameters<T> &>(ach);
-	BFGS_hyperparameters<T> *bach = dynamic_cast<BFGS_hyperparameters<T> *>(achieved);
-	if(bach) *bach = ach;
-	}
-
-return optresult<T>(f,converged);
+return st.finish(achieved);
 }
 
 
 
 //force instantiation
-template double linesearch_wolfe(const OptProblem<double> *o, NRVec<double> &x, double &f, NRVec<double> &g, const NRVec<double> &p, double alpha1, const double c1, const double c2, const double alphamax, const int maxiter, int &nfunc, const int verbose);
+template double linesearch_wolfe(optevaluator<double> &e, NRVec<double> &x, double &f, NRVec<double> &g, const NRVec<double> &p, const double alpha1, const linesearch_hyperparameters<double> &lp);
+template class optstate<double>;
 template class testproblem<double>;
 template class BFGS_Optimizer<double>;
 
diff --git a/optimize.h b/optimize.h
index 64439ab..da3fbef 100644
--- a/optimize.h
+++ b/optimize.h
@@ -52,8 +52,12 @@ virtual ~testproblem() {};
 
 //base classes for optimization methods
 
+template<typename T> class optstate;
+
+
 //this is a collection of parameters, which can be augmented by method-specific additions, to be used as a base class
-//methods obtain their specific parameters via dynamic_cast, missing ones are taken from the default constructor
+//methods take their specific parameters via setfrom(), those not provided by the caller keep the defaults from the constructor
+//each derived class has to override setfrom(), converged() and setreached() if it adds parameters relevant for them
 template<typename T>
 class hyperparameters {
 public:
@@ -61,8 +65,12 @@ int my_maxstep; //maximum number of iterations
 int my_maxfunc; //maximum number of function evaluations, <=0 means unlimited
 T my_f_thr; //threshhold on function change in a single step
 T my_x_thr; //threshold on x norm change in a single step
+int my_saddle_order; //0=minimum, k=saddle point with k negative Hessian eigenvalues (only for methods supporting it)
 int my_verbose; //0=silent, 1=per-iteration info, 2=more detailed info
-hyperparameters() : my_maxstep(200), my_maxfunc(0), my_f_thr(0), my_x_thr(0), my_verbose(1) {};
+hyperparameters() : my_maxstep(200), my_maxfunc(0), my_f_thr(0), my_x_thr(0), my_saddle_order(0), my_verbose(1) {};
+virtual void setfrom(const hyperparameters &p) {*this=p;} //copy those parameters from p which exist in both classes
+virtual bool converged(const optstate<T> &s) const {return s.stepdone && s.df<my_f_thr && s.dx<my_x_thr;}
+virtual void setreached(const optstate<T> &s) {my_maxstep=s.iter; my_maxfunc=s.e.nfunc; my_f_thr=s.df; my_x_thr=s.dx;} //store the values reached in the minimization
 virtual ~hyperparameters(){};
 };
 
@@ -72,6 +80,21 @@ class gradient_hyperparameters : public hyperparameters<T> {
 public:
 T my_g_thr; //gradient norm convergence threshold
 gradient_hyperparameters() : my_g_thr(1e-6) {};
+virtual void setfrom(const hyperparameters<T> &p) {const gradient_hyperparameters *q=dynamic_cast<const gradient_hyperparameters *>(&p); if(q) *this= *q; else hyperparameters<T>::setfrom(p);}
+virtual bool converged(const optstate<T> &s) const {return (s.gnorm>=0 && s.gnorm<my_g_thr) || hyperparameters<T>::converged(s);}
+virtual void setreached(const optstate<T> &s) {hyperparameters<T>::setreached(s); my_g_thr=s.gnorm;}
+};
+
+//parameters of gradient methods using line search
+template<typename T>
+class linesearch_hyperparameters : public gradient_hyperparameters<T> {
+public:
+T my_c1; //sufficient decrease (Armijo) parameter of Wolfe conditions, 0<c1<c2<1
+T my_c2; //curvature parameter of Wolfe conditions
+int my_maxlinesearch; //max. number of function evaluations in one line search
+T my_maxsteplength; //max. norm of a single step in x
+linesearch_hyperparameters() : my_c1(1e-4), my_c2(0.9), my_maxlinesearch(30), my_maxsteplength(1e10) {};
+virtual void setfrom(const hyperparameters<T> &p) {const linesearch_hyperparameters *q=dynamic_cast<const linesearch_hyperparameters *>(&p); if(q) *this= *q; else gradient_hyperparameters<T>::setfrom(p);}
 };
 
 
@@ -86,6 +109,48 @@ operator T() const {return f;}
 };
 
 
+//evaluation of the optimized function, counts the evaluations and keeps the best point found
+template<typename T>
+class optevaluator {
+public:
+const OptProblem<T> *o;
+int maxfunc; //<=0 means unlimited
+int nfunc; //no. of evaluations done
+T fbest; //lowest function value found
+NRVec<T> xbest; //and its point
+optevaluator(const OptProblem<T> *oo, const int maxf=0) : o(oo), maxfunc(maxf), nfunc(0), fbest(0) {};
+T operator()(const NRVec<T> &x, NRVec<T> *grad=NULL)
+	{
+	T f=o->funcgrad(x,grad);
+	if(nfunc++ ==0 || f<fbest) {fbest=f; xbest=x;}
+	return f;
+	}
+bool exhausted() const {return maxfunc>0 && nfunc>=maxfunc;}
+int remaining() const {return maxfunc>0 ? maxfunc-nfunc : std::numeric_limits<int>::max();}
+};
+
+
+//bookkeeping of the iterations common to all methods: convergence and limit tests, verbose output, results
+template<typename T>
+class optstate {
+public:
+const char *name; //of the method for verbose output
+const hyperparameters<T> &hp; //parameters used
+const optevaluator<T> &e;
+int iter; //no. of iterations done
+T f; //current function value
+T gnorm; //current gradient norm, negative if not available
+T df; //|f change| in the last step
+T dx; //|x change| in the last step
+bool stepdone; //df and dx are valid
+bool converged;
+optstate(const char *nm, const hyperparameters<T> &h, const optevaluator<T> &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
+void step(const T dff, const T dxx) {df=std::abs(dff); dx=dxx; stepdone=true;} //record a step done
+optresult<T> finish(hyperparameters<T> *achieved) const; //final report, fills achieved if not NULL
+};
+
+
 //base class for optimization methods, virtual functions to allow a pointer to derived class
 template<typename T>
 class Optimizer {
@@ -102,22 +167,19 @@ virtual ~Optimizer() {};
 
 //line search along direction p satisfying strong Wolfe conditions (cf. Nocedal & Wright, Numerical Optimization, 2nd ed., alg. 3.5 and 3.6)
 //on input x,f,g at the starting point, on successful output x,f,g at the new point; returns step length or 0 on failure (x,f,g then unchanged)
-//maxiter = max. number of function evaluations, nfunc is incremented by the number of evaluations done
+//alpha1 = initial trial step length
 template<typename T>
-T linesearch_wolfe(const OptProblem<T> *o, NRVec<T> &x, T &f, NRVec<T> &g, const NRVec<T> &p, T alpha1, const T c1, const T c2, const T alphamax, const int maxiter, int &nfunc, const int verbose);
+T linesearch_wolfe(optevaluator<T> &e, NRVec<T> &x, T &f, NRVec<T> &g, const NRVec<T> &p, const T alpha1, const linesearch_hyperparameters<T> &lp);
 
 
 //specific methods
 
 template<typename T>
-class BFGS_hyperparameters : public gradient_hyperparameters<T> {
+class BFGS_hyperparameters : public linesearch_hyperparameters<T> {
 public:
-T my_c1; //sufficient decrease (Armijo) parameter of Wolfe conditions, 0<c1<c2<1
-T my_c2; //curvature parameter of Wolfe conditions
-int my_maxlinesearch; //max. number of function evaluations in one line search
-T my_maxsteplength; //max. norm of a single step in x
 T my_sy_thr; //skip the update if s.y <= sy_thr*|s|*|y| (curvature condition violated numerically)
-BFGS_hyperparameters() : my_c1(1e-4), my_c2(0.9), my_maxlinesearch(30), my_maxsteplength(1e10), my_sy_thr(1e-12) {};
+BFGS_hyperparameters() : my_sy_thr(1e-12) {};
+virtual void setfrom(const hyperparameters<T> &p) {const BFGS_hyperparameters *q=dynamic_cast<const BFGS_hyperparameters *>(&p); if(q) *this= *q; else linesearch_hyperparameters<T>::setfrom(p);}
 };
 
 template<typename T>