float support for optimize

optimize.cc

@@ -1124,7 +1124,7 @@ NRSMat<T> A=J.transposedtimes(); //J^T J
 NRVec<T> g=r*J; //J^T r, the gradient of f
 NRVec<T> diag(n);
 for(int i=0; i<n; ++i) diag[i]=((const NRSMat<T> &)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<T>::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<T> hh(h); hh*=h;
+		NRVec<T> xx(x); xx*=x; xx+=eps*eps;
+		hh/=xx;
+		if(hh.amax()<=100*eps*eps)
 			{
 			if(verbose) std::cout <<name<<": step at machine precision\n";
 			st.converged=true;
@@ -2080,6 +2085,7 @@ template optresult<T> firstorder_minimize(const char *name, firstorder_step<T> &
 template optresult<T> trustregion_minimize(const char *name, trustregion_step<T> &stepper, const OptProblem<T> *o, NRVec<T> &x, const trustregion_hyperparameters<T> &hp, hyperparameters<T> *achieved); \
 
 INSTANTIZE(double)
+INSTANTIZE(float)
 
 
 }//namespace

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<typename T> struct opt_precision;
+template<> struct opt_precision<double>
+	{
+	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<float>
+	{
+	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<typename T>
 class OptLeastSquares : public OptProblem<T>, public OptResiduals<T> {
 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<T>::fd_forward_step()) : fd_step(h) {};
 void jacobian_fd(const NRVec<T> &x, const NRVec<T> &r, NRMat<T> &J) const; //Jacobian by forward differences, r = residuals at x
 virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const;
 virtual ~OptLeastSquares() {};
@@ -164,7 +195,7 @@ class ParameterShift : public OptProblemAdapter<T> {
 public:
 const OptTrigonometric<T> *trig;
 T h;
-ParameterShift(const OptProblem<T> *pp, const T hh=1e-6);
+ParameterShift(const OptProblem<T> *pp, const T hh=opt_precision<T>::fd_central_step());
 virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const;
 };
 
@@ -173,7 +204,7 @@ template<typename T>
 class FiniteDifference : public OptProblemAdapter<T> {
 public:
 T h;
-FiniteDifference(const OptProblem<T> *pp, const T hh=1e-6) : OptProblemAdapter<T>(pp), h(hh) {};
+FiniteDifference(const OptProblem<T> *pp, const T hh=opt_precision<T>::fd_central_step()) : OptProblemAdapter<T>(pp), h(hh) {};
 virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const;
 };
 
@@ -207,7 +238,7 @@ template<typename T>
 class gradient_hyperparameters : public hyperparameters<T> {
 public:
 T my_g_thr; //gradient norm convergence threshold
-gradient_hyperparameters() : my_g_thr(1e-6) {};
+gradient_hyperparameters() : my_g_thr(opt_precision<T>::g_thr()) {};
 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;}
@@ -379,7 +410,7 @@ template<typename T>
 class quasinewton_hyperparameters : public linesearch_hyperparameters<T> {
 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<T>::sy_thr()) {};
 virtual void setfrom(const hyperparameters<T> &p) {const quasinewton_hyperparameters *q=dynamic_cast<const quasinewton_hyperparameters *>(&p); if(q) *this= *q; else linesearch_hyperparameters<T>::setfrom(p);}
 };
 
@@ -392,7 +423,7 @@ class BFGS_hyperparameters : public quasinewton_hyperparameters<T> {
 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<T>::sr1_thr()) {};
 virtual void setfrom(const hyperparameters<T> &p) {const BFGS_hyperparameters *q=dynamic_cast<const BFGS_hyperparameters *>(&p); if(q) *this= *q; else quasinewton_hyperparameters<T>::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<T>::update_thr()) {};
 virtual void setfrom(const hyperparameters<T> &p) {const trustregion_hyperparameters *q=dynamic_cast<const trustregion_hyperparameters *>(&p); if(q) *this= *q; else gradient_hyperparameters<T>::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<T>::fd_forward_step()) {};
 virtual void setfrom(const hyperparameters<T> &p) {const LM_hyperparameters *q=dynamic_cast<const LM_hyperparameters *>(&p); if(q) *this= *q; else gradient_hyperparameters<T>::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<T>::nm_f_thr(); this->my_x_thr=opt_precision<T>::nm_x_thr();};
 virtual void setfrom(const hyperparameters<T> &p) {const NelderMead_hyperparameters *q=dynamic_cast<const NelderMead_hyperparameters *>(&p); if(q) *this= *q; else hyperparameters<T>::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<T>::powell_f_thr(); this->my_x_thr=opt_precision<T>::powell_x_thr();};
 virtual void setfrom(const hyperparameters<T> &p) {const Powell_hyperparameters *q=dynamic_cast<const Powell_hyperparameters *>(&p); if(q) *this= *q; else hyperparameters<T>::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<T>::rhoend()), my_flexible_radius(true) {this->my_maxstep=100000;};
 virtual void setfrom(const hyperparameters<T> &p) {const Cobyla_hyperparameters *q=dynamic_cast<const Cobyla_hyperparameters *>(&p); if(q) *this= *q; else hyperparameters<T>::setfrom(p);}
 };
 

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<float>) on float test problems
+template<class O, class H>
+static float runf(O &o, const OptProblem<float> *p, NRVec<float> &x, H &h, bool *converged)
+{
+h.my_verbose=0;
+x.copyonwrite();
+optresult<float> r=o.minimize(p,x,&h,NULL);
+*converged=r.converged;
+return r.f;
+}
+
+static NRVec<float> startpointf(const int n, const int seed)
+{
+NRVec<double> xd=startpoint(n,seed);
+NRVec<float> x(n);
+for(int i=0; i<n; ++i) x[i]=xd[i];
+return x;
+}
+
+//the same method in double and in float from the same start, the float run with its default thresholds
+template<template<class> class O, template<class> class H>
+static void compare_float(const char *name, const OptProblem<double> *pd, const OptProblem<float> *pf, const int n, const int seed, const int maxstep, void (*modify)(H<double> &, H<float> &)=NULL)
+{
+NRVec<double> xd=startpoint(n,seed); NRVec<float> xf(n); for(int i=0; i<n; ++i) xf[i]=xd[i];
+H<double> hd; H<float> hf; hd.my_maxstep=hf.my_maxstep=maxstep; if(modify) modify(hd,hf);
+O<double> od; O<float> of;
+bool convd, convf;
+hd.my_verbose=0; xd.copyonwrite(); optresult<double> 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<double> &hd, BFGS_hyperparameters<float> &hf) {hd.my_update=hf.my_update=QN_DFP; hd.my_c2=hf.my_c2=0.1;}
+static void sr1(BFGS_hyperparameters<double> &hd, BFGS_hyperparameters<float> &hf) {hd.my_update=hf.my_update=QN_SR1;}
+
+static void test_float()
+{
+cout <<"--- single precision (float)\n";
+testproblem<float> 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<double> vd; testproblem_vqe<float> vf;
+compare_float<BFGS_Optimizer,BFGS_hyperparameters>("BFGS on testproblem_vqe",&vd,&vf,8,3,5000);
+compare_float<BFGS_Optimizer,BFGS_hyperparameters>("DFP on testproblem_vqe",&vd,&vf,8,3,5000,&dfp);
+compare_float<BFGS_Optimizer,BFGS_hyperparameters>("SR1 on testproblem_vqe",&vd,&vf,8,3,5000,&sr1);
+compare_float<LBFGS_Optimizer,LBFGS_hyperparameters>("L-BFGS on testproblem_vqe",&vd,&vf,8,3,5000);
+compare_float<CG_Optimizer,CG_hyperparameters>("CG on testproblem_vqe",&vd,&vf,8,3,5000);
+compare_float<GDIIS_Optimizer,GDIIS_hyperparameters>("GDIIS on testproblem_vqe",&vd,&vf,8,3,5000);
+{NelderMead_hyperparameters<float> h; h.my_maxstep=100000; NelderMead_Optimizer<float> o; NRVec<float> 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<float> h; Powell_Optimizer<float> o; NRVec<float> 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<float> h; Cobyla_Optimizer<float> o; NRVec<float> 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<float> h; h.my_maxstep=20000; h.my_lr=1e-2; Adam_Optimizer<float> o; NRVec<float> 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<float> fd(&tp); NRVec<float> 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<double> fdd(&vd); FiniteDifference<float> fdf(&vf);
+compare_float<BFGS_Optimizer,BFGS_hyperparameters>("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<double> xd=startpoint(6,3); NRVec<float> xf(6); for(int i=0;i<6;++i) xf[i]=xd[i];
+TrustRegion_hyperparameters<double> hd; TrustRegion_Optimizer<double> od; double fd=run(od,&vd,xd,hd);
+TrustRegion_hyperparameters<float> hf; TrustRegion_Optimizer<float> 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<float> p; NRVec<float> b(2); for(int i=0;i<2;++i) b[i]=p.start(0,i);
+LM_hyperparameters<float> h; h.my_maxstep=100000; h.my_scaling=true; LM_Optimizer<float> 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 <<"=== "<<npass<<" passed, "<<nfail<<" failed\n";
 return nfail;
 }