From b38bd343f464ae8a8fa6b560a6e83848dbfd10a2 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Mon, 26 Oct 2026 16:10:07 +0100 Subject: [PATCH] testoptimize testbed created --- .gitignore | 1 + Makefile.am | 6 +- optimize.cc | 118 -------- optimize.h | 34 --- t.cc | 1 + testoptimize.cc | 733 ++++++++++++++++++++++++++++++++++++++++++++++++ testoptimize.h | 178 ++++++++++++ 7 files changed, 917 insertions(+), 154 deletions(-) create mode 100644 testoptimize.cc create mode 100644 testoptimize.h diff --git a/.gitignore b/.gitignore index 55eccbd..3613e92 100644 --- a/.gitignore +++ b/.gitignore @@ -65,4 +65,5 @@ gmon.out npytest *.npy *.prompt +testoptimize # CVS default ignores end diff --git a/Makefile.am b/Makefile.am index 396a8f6..700e944 100644 --- a/Makefile.am +++ b/Makefile.am @@ -2,13 +2,15 @@ lib_LTLIBRARIES = libla.la include_HEADERS = LA_version.h optimize.h tensor.h miscfunc.h simple.h vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h davidson.h lanczos.h laerror.h mat.h qsort.h vec.h bisection.h diis.h la.h noncblas.h smat.h numbers.h bitvector.h fourindex.h la_traits.h la_random.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h reg.h regsurf.h libla_la_SOURCES = simple.cc optimize.cc tensor.cc miscfunc.cc quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc numbers.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc contfrac.cc graph.cc la_random.cc reg.cc regsurf.cc nodist_libla_la_SOURCES = version.cc -check_PROGRAMS = t test tX test_reg test_regsurf -t_SOURCES = t.cc t2.cc +noinst_HEADERS = testoptimize.h +check_PROGRAMS = t test tX test_reg test_regsurf testoptimize +t_SOURCES = t.cc t2.cc #nodist_t_SOURCES = test_SOURCES = test.cc tX_SOURCES = tX.cc test_reg_SOURCES = test_reg.cc test_regsurf_SOURCES = test_regsurf.cc +testoptimize_SOURCES = testoptimize.cc LDADD = .libs/libla.a ACLOCAL_AMFLAGS = -I m4 EXTRA_DIST = get_git_version LICENSE doxygen.cfg aminclude.am acinclude.m4 footer.html diff --git a/optimize.cc b/optimize.cc index aa3f599..f17ce6a 100644 --- a/optimize.cc +++ b/optimize.cc @@ -24,40 +24,6 @@ namespace LA { -template -T testproblem::funcgrad(const NRVec &x, NRVec *grad) const -{ -int n=x.size(); -T r= n*(n-1)/2; -for(int i=0; isize()!=x.size()) laerror("grad size mismatch"); - grad->copyonwrite(); //I forgot this initially - for(int k=0; k -void testproblem_vqe_constrained::constraints(const NRVec &x, NRVec &c) const -{ -c.copyonwrite(); -c[0]=0; -for(int i=0; i static T central_difference(const OptProblem *p, NRVec &xs, const NRVec &x, const int i, const T h) @@ -120,87 +86,6 @@ return f; } -//Bloch vector components z = cos a cos b, x = cos a sin b of the qubit RY(b) RX(a)|0> and their derivatives -template -class blochqubit { -public: -T z, x, za, zb, xa, xb; -blochqubit(const T a, const T b) - { - T ca=std::cos(a), sa=std::sin(a), cb=std::cos(b), sb=std::sin(b); - z=ca*cb; x=ca*sb; za= -sa*cb; zb= -ca*sb; xa= -sa*sb; xb=ca*cb; - } -}; - -template -T testproblem_vqe::funcgrad(const NRVec &x, NRVec *grad) const -{ -const int nq=x.size()/2; -if(2*nq!=x.size()) laerror("testproblem_vqe needs an even number of parameters"); -T r=0; -for(int i=0; i q(x[2*i],x[2*i+1]); - r+= hfield(i)*q.z + kfield(i)*q.x; - if(isize()!=x.size()) laerror("grad size mismatch"); - grad->copyonwrite(); - for(int i=0; i q(x[2*i],x[2*i+1]); - T eff=hfield(i); //coefficient of z_i including the couplings - if(i>0) eff+= J(i-1)*std::cos(x[2*i-2])*std::cos(x[2*i-1]); - if(i -void testproblem_vqe::hessian(const NRVec &x, NRSMat &h) const -{ -const int nq=x.size()/2; -h=(T)0; -for(int i=0; i q(x[2*i],x[2*i+1]); - T eff=hfield(i); - if(i>0) eff+= J(i-1)*std::cos(x[2*i-2])*std::cos(x[2*i-1]); - if(i q1(x[2*i+2],x[2*i+3]); - h(2*i,2*i+2)= J(i)*q.za*q1.za; - h(2*i,2*i+3)= J(i)*q.za*q1.zb; - h(2*i+1,2*i+2)= J(i)*q.zb*q1.za; - h(2*i+1,2*i+3)= J(i)*q.zb*q1.zb; - } - } -} - -template -void testproblem_vqe::metric(const NRVec &x, NRSMat &g) const -{ -const int nq=x.size()/2; -g=(T)0; -for(int i=0; i static T cubic_interpolate(const T a, const T fa, const T da, const T b, const T fb, const T db) @@ -1666,9 +1551,6 @@ 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; -template class testproblem; -template class testproblem_vqe; -template class testproblem_vqe_constrained; template class ParameterShift; template class FiniteDifference; template class BFGS_Optimizer; diff --git a/optimize.h b/optimize.h index 33fbf0a..ef85e8c 100644 --- a/optimize.h +++ b/optimize.h @@ -125,40 +125,6 @@ virtual T funcgrad(const NRVec &x, NRVec *grad) const; }; -//test case -template -class testproblem : public OptProblem { -public: -virtual T funcgrad(const NRVec &x, NRVec *grad) const; -virtual ~testproblem() {}; -}; - -//test case with all capabilities: Ising chain with longitudinal and transverse fields in the product state of qubits RY(b_i) RX(a_i)|0>, -//x = (a_0,b_0,a_1,b_1,...), E = sum_i (h_i z_i + k_i x_i) + sum_i J_i z_i z_{i+1}, Bloch vector components z_i = cos a_i cos b_i, x_i = cos a_i sin b_i; -//trigonometric of degree 1 with period 2 pi in all parameters, Fubini-Study metric 1/4 diag(1, cos^2 a_i) for each qubit -template -class testproblem_vqe : public OptProblem, public OptHessian, public OptMetric, public OptTrigonometric { -public: -static T hfield(const int i) {return 0.7+0.1*(i%3);} -static T kfield(const int i) {return 0.4;} -static T J(const int i) {return 1;} -virtual T funcgrad(const NRVec &x, NRVec *grad) const; -virtual void hessian(const NRVec &x, NRSMat &h) const; -virtual void metric(const NRVec &x, NRSMat &g) const; -virtual int degree(const int i) const {return 1;} -virtual ~testproblem_vqe() {}; -}; - -//constrained test case: testproblem_vqe with nonnegative total magnetization sum_i z_i >= 0 -template -class testproblem_vqe_constrained : public testproblem_vqe, public OptConstraints { -public: -virtual int nconstraints() const {return 1;} -virtual void constraints(const NRVec &x, NRVec &c) const; -virtual ~testproblem_vqe_constrained() {}; -}; - - //base classes for optimization methods template class optstate; diff --git a/t.cc b/t.cc index de77c6a..f9311aa 100644 --- a/t.cc +++ b/t.cc @@ -20,6 +20,7 @@ #include #include "la.h" +#include "testoptimize.h" #include "vecmat3.h" #include "quaternion.h" #include "permutation.h" diff --git a/testoptimize.cc b/testoptimize.cc new file mode 100644 index 0000000..e44a0e7 --- /dev/null +++ b/testoptimize.cc @@ -0,0 +1,733 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2026 Jiri Pittner or + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +//regression tests of the optimization methods and of the functions of the vector and matrix classes they rely on +//each check prints PASS or FAIL, the exit code is the number of failures + +#include +#include +#include +#include +#include "la.h" +#include "testoptimize.h" +#include "diis.h" + +using namespace std; +using namespace LA; + +static int npass=0, nfail=0; + +static void check(const char *what, const bool ok, const double value=0) +{ +if(ok) ++npass; else ++nfail; +cout <<(ok?"PASS ":"FAIL ")< { +public: +double funcgrad(const NRVec &x, NRVec *g) const + { + int n=x.size(); + double r=0; + if(g) {g->copyonwrite(); *g=0.;} + for(int i=0; i { +public: +double funcgrad(const NRVec &x, NRVec *g) const + { + int n=x.size(); + double r=0; + if(g) {g->copyonwrite(); *g=0.;} + for(int i=0; i { +public: +double sigma; +noisyquadratic(const double s) : sigma(s) {}; +static double exact(const NRVec &x) {double r=0; for(int i=0; i &x, NRVec *g) const + { + if(g) {g->copyonwrite(); for(int i=0; i, public OptTrigonometric { +public: +static double tr(const double t) {double w=2*M_PI/5; return 0.3+cos(w*t)-0.7*sin(w*t)+0.4*cos(2*w*t)+0.2*sin(2*w*t)-0.5*cos(3*w*t)+0.9*sin(3*w*t);} +static double dtr(const double t) {double w=2*M_PI/5; return w*(-sin(w*t)-0.7*cos(w*t)-0.8*sin(2*w*t)+0.4*cos(2*w*t)+1.5*sin(3*w*t)+2.7*cos(3*w*t));} +double funcgrad(const NRVec &x, NRVec *g) const + { + double s=1.5+0.5*tanh(x[1]); + if(g) {g->copyonwrite(); (*g)[0]=dtr(x[0])*s; double th=tanh(x[1]); (*g)[1]=tr(x[0])*0.5*(1-th*th)+2*(x[1]-0.3);} + return tr(x[0])*s+(x[1]-0.3)*(x[1]-0.3); + } +int degree(const int i) const {return i==0 ? 3 : 0;} +double period(const int i) const {return 5;} +}; + +//counts the evaluations of the wrapped problem, its capabilities remain accessible +class counting : public OptProblemAdapter { +public: +mutable int count; +counting(const OptProblem *pp) : OptProblemAdapter(pp), count(0) {}; +double funcgrad(const NRVec &x, NRVec *g) const {++count; return this->p->funcgrad(x,g);} +}; + +//constrained classics with known solutions +class quadratic_linearconstraint : public OptProblem, public OptConstraints { //min (x0-1)^2+(x1-2)^2, x0+x1<=1: f=2 at (0,1) +public: +double funcgrad(const NRVec &x, NRVec *g) const {return (x[0]-1)*(x[0]-1)+(x[1]-2)*(x[1]-2);} +int nconstraints() const {return 1;} +void constraints(const NRVec &x, NRVec &c) const {c.copyonwrite(); c[0]=1-x[0]-x[1];} +}; +class product_disk : public OptProblem, public OptConstraints { //min x0 x1, x0^2+x1^2<=1: f=-1/2 (Powell 1994) +public: +double funcgrad(const NRVec &x, NRVec *g) const {return x[0]*x[1];} +int nconstraints() const {return 1;} +void constraints(const NRVec &x, NRVec &c) const {c.copyonwrite(); c[0]=1-x[0]*x[0]-x[1]*x[1];} +}; +class rosensuzuki : public OptProblem, public OptConstraints { //f=-44 at (0,1,2,-1) +public: +double funcgrad(const NRVec &x, NRVec *g) const {return x[0]*x[0]+x[1]*x[1]+2*x[2]*x[2]+x[3]*x[3]-5*x[0]-5*x[1]-21*x[2]+7*x[3];} +int nconstraints() const {return 3;} +void constraints(const NRVec &x, NRVec &c) const + { + c.copyonwrite(); + c[0]=8-x[0]*x[0]-x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-x[0]+x[1]-x[2]+x[3]; + c[1]=10-x[0]*x[0]-2*x[1]*x[1]-x[2]*x[2]-2*x[3]*x[3]+x[0]+x[3]; + c[2]=5-2*x[0]*x[0]-x[1]*x[1]-x[2]*x[2]-2*x[0]+x[1]+x[3]; + } +}; +class sum_circle : public OptProblem, public OptConstraints { //min x0+x1, x0^2+x1^2=1 as two inequalities: f=-sqrt(2) +public: +double funcgrad(const NRVec &x, NRVec *g) const {return x[0]+x[1];} +int nconstraints() const {return 2;} +void constraints(const NRVec &x, NRVec &c) const {c.copyonwrite(); c[0]=x[0]*x[0]+x[1]*x[1]-1; c[1]= -c[0];} +}; + + +//auxiliary functions + +static NRVec startpoint(const int n, const int seed, const double scale=1) +{ +NRVec x(n); +srand(seed); +x.randomize(scale); +return x; +} + +static NRVec rosenbrockstart(const int n) +{ +NRVec x(n); +for(int i=0; i &p, const NRVec &x) +{ +const int n=x.size(); +NRVec g(n); +p.funcgrad(x,&g); +double err=0; +for(int i=0; i a(x), b(x); + a.copyonwrite(); b.copyonwrite(); + a[i]+=1e-6; b[i]-=1e-6; + err=max(err,abs((p.funcgrad(a,NULL)-p.funcgrad(b,NULL))/2e-6-g[i])); + } +return err; +} + +//minimizes and returns f, x is modified, achieved filled +static double run(Optimizer &o, const OptProblem *p, NRVec &x, hyperparameters &h, hyperparameters *achieved=NULL, bool *converged=NULL) +{ +h.my_verbose=0; +x.copyonwrite(); +optresult r=o.minimize(p,x,&h,achieved); +if(converged) *converged=r.converged; +return r.f; +} + +//number of maxfunc limits from 1 to maxmax for which the evaluation budget is exceeded or misreported (from minfunc on, below it the initial evaluations are unavoidable) +static int maxfunc_violations(Optimizer &o, const OptProblem *p, const NRVec &x0, hyperparameters &h, const int minfunc, const int maxmax) +{ +int bad=0; +for(int mf=minfunc; mf<=maxmax; ++mf) + { + counting cp(p); + NRVec x(x0); + h.my_maxfunc=mf; + hyperparameters a; + run(o,&cp,x,h,&a); + if(cp.count>mf || a.my_maxfunc!=cp.count) ++bad; + } +return bad; +} + + +//numerical Fubini-Study metric of one qubit RY(b) RX(a)|0> +static void fubinistudy(const double a, const double b, double g[2][2]) +{ +typedef complex C; +const double h=1e-5; +C psi[2], d[2][2]; +for(int k=-1; k<=1; ++k) //k=-1: the state, k=0,1: derivatives + { + C p[2], m[2]; + for(int s= -1; s<=1; s+=2) + { + double aa= a+(k==0 ? s*h : 0), bb= b+(k==1 ? s*h : 0); + C c0=cos(aa/2), c1=C(0,-sin(aa/2)); + C *q= s>0 ? p : m; + q[0]=cos(bb/2)*c0-sin(bb/2)*c1; + q[1]=sin(bb/2)*c0+cos(bb/2)*c1; + } + if(k<0) {psi[0]=(p[0]+m[0])/2.; psi[1]=(p[1]+m[1])/2.;} + else for(int s=0; s<2; ++s) d[k][s]=(p[s]-m[s])/(2*h); + } +for(int i=0; i<2; ++i) for(int j=0; j<2; ++j) + { + C dd=conj(d[i][0])*d[j][0]+conj(d[i][1])*d[j][1]; + C di=conj(d[i][0])*psi[0]+conj(d[i][1])*psi[1]; + C dj=conj(psi[0])*d[j][0]+conj(psi[1])*d[j][1]; + g[i][j]=real(dd)-real(di*dj); + } +} + + +//reference implementation of the PyTorch update formulas of Adam, AdamW, AMSGrad and SGD with momentum, for comparison +static NRVec reference_firstorder(const OptProblem &p, NRVec x, const int steps, const double lr, const int type, const double wd) +{ +//type 0 Adam/AdamW, 1 AMSGrad, 2 SGD, 3 SGD-Nesterov +const int n=x.size(); +const double b1=0.9, b2=0.999, eps=1e-8, mom=0.9; +x.copyonwrite(); +NRVec g(n), m(0.,n), v(0.,n), vmax(0.,n), vel(0.,n); +p.funcgrad(x,&g); +for(int k=1; k<=steps; ++k) + { + for(int i=0; i gn(n); + p.funcgrad(x,&gn); + g=gn; + } +return x; +} + + +static void test_problems_and_capabilities() +{ +cout <<"--- test problems, capabilities and gradient adapters\n"; +testproblem tp; +testproblem_vqe vqe; +NRVec x=startpoint(8,4,2.); +check("testproblem gradient vs finite differences",fd_gradient_error(tp,x)<1e-8,fd_gradient_error(tp,x)); +check("testproblem_vqe gradient vs finite differences",fd_gradient_error(vqe,x)<1e-8,fd_gradient_error(vqe,x)); +{ +NRSMat h(8), hfd(8); +vqe.hessian(x,h); +for(int i=0; i<8; ++i) + { + NRVec a(x), b(x), ga(8), gb(8); + a.copyonwrite(); b.copyonwrite(); + a[i]+=1e-6; b[i]-=1e-6; + vqe.funcgrad(a,&ga); vqe.funcgrad(b,&gb); + for(int j=0; j<=i; ++j) hfd(i,j)=(ga[j]-gb[j])/2e-6; + } +NRSMat d=h-hfd; +check("testproblem_vqe Hessian vs finite differences of the gradient",d.norm()<1e-8,d.norm()); +NRVec v=startpoint(8,5), hv=vqe.hessian_times(x,v), hv2=h*v; +hv-=hv2; +check("OptHessian::hessian_times vs hessian*v",hv.norm()<1e-14,hv.norm()); +} +{ +NRSMat g(8); +vqe.metric(x,g); +double err=0; +for(int q=0; q<4; ++q) + { + double fs[2][2]; + fubinistudy(x[2*q],x[2*q+1],fs); + for(int i=0; i<2; ++i) for(int j=0; j<2; ++j) err=max(err,abs(g(2*q+i,2*q+j)-fs[i][j])); + } +check("testproblem_vqe metric vs numerical Fubini-Study metric of the state",err<1e-9,err); +} +{ +NRVec g(8), gps(8), gfd(8); +vqe.funcgrad(x,&g); +ParameterShift ps(&vqe); +FiniteDifference fd(&vqe); +ps.funcgrad(x,&gps); fd.funcgrad(x,&gfd); +gps-=g; gfd-=g; +check("ParameterShift gradient, degree 1",gps.norm()<1e-13,gps.norm()); +check("FiniteDifference gradient",gfd.norm()<1e-8,gfd.norm()); +trigdegree3 t3; +ParameterShift ps3(&t3); +NRVec y(2), g3(2), g3ps(2); +y[0]=1.234; y[1]=0.6; +t3.funcgrad(y,&g3); ps3.funcgrad(y,&g3ps); +g3ps-=g3; +check("ParameterShift gradient, degree 3 and period 5, and a non-trigonometric parameter",g3ps.norm()<1e-8,g3ps.norm()); +FiniteDifference fdps(&ps); +check("optcapability through adapters",optcapability >(&ps)==&vqe && optcapability >(&fdps)==&vqe && optcapability >(&tp)==NULL); +testproblem_vqe_constrained vc; +counting cvc(&vc); +check("constraints found through an adapter",optcapability >(&cvc)==&vc); +} +} + + +static void test_linemin() +{ +cout <<"--- derivative-free line minimization (Brent)\n"; +class oned : public OptProblem {public: int k; double funcgrad(const NRVec &x, NRVec *g) const {double t=x[0]; + switch(k) {case 0: return (t-3.7)*(t-3.7)+1; case 1: return cos(t); case 2: return exp(t)-5*t; case 3: return pow(t+2.5,4); default: return abs(t-0.3);}}}; +const double tmin[]={3.7,M_PI,log(5.),-2.5,0.3}; +bool ok=true; +double worst=0; +for(int k=0; k<5; ++k) for(double t1: {0.1,1.}) + { + oned p; p.k=k; + optevaluator e(&p); + NRVec x(0.,1), u(1.,1); + double f=e(x); + linemin_brent(e,x,f,u,t1,1e-8,1e-10,100); + double err=abs(x[0]-tmin[k]); + if(err>(k==3 ? 1e-4 : 1e-6)) ok=false; + worst=max(worst,err); + } +check("linemin_brent finds the minima of 5 one-dimensional functions",ok,worst); +oned p; p.k=0; +bool budget=true; +for(int mb=1; mb<=10; ++mb) {optevaluator e(&p); NRVec x(0.,1), u(1.,1); double f=e(x); linemin_brent(e,x,f,u,0.1,1e-8,1e-10,mb); if(e.nfunc-1>mb) budget=false;} +check("linemin_brent respects the evaluation budget",budget); +} + + +static void test_gradient_methods() +{ +cout <<"--- gradient methods\n"; +testproblem tp; +rosenbrock rb; +//BFGS, DFP, SR1 +const QN_update upd[]={QN_BFGS,QN_DFP,QN_SR1}; +const char *updname[]={"BFGS","DFP","SR1"}; +for(int u=0; u<3; ++u) + { + BFGS_hyperparameters h; h.my_update=upd[u]; h.my_g_thr=1e-6; h.my_maxstep=5000; + if(upd[u]==QN_DFP) h.my_c2=0.1; //DFP needs an accurate line search + BFGS_Optimizer o; + NRVec x=startpoint(10,1); + bool conv; + double f=run(o,&tp,x,h,NULL,&conv); + string s=string(updname[u])+" converges on testproblem n=10"; + check(s.c_str(),conv && f<1e-10,f); + x=rosenbrockstart(10); + f=run(o,&rb,x,h,NULL,&conv); + s=string(updname[u])+" converges on Rosenbrock n=10"; + check(s.c_str(),conv && f<1e-10,f); + } +{ +LBFGS_hyperparameters h; h.my_g_thr=1e-6; h.my_maxstep=10000; +LBFGS_Optimizer o; +NRVec x=rosenbrockstart(100); +bool conv; +double f=run(o,&rb,x,h,NULL,&conv); +check("L-BFGS converges on Rosenbrock n=100",conv && f<1e-10,f); +} +const CG_beta beta[]={CG_FR,CG_PR,CG_HS,CG_DY,CG_HZ}; +for(int b=0; b<5; ++b) + { + CG_hyperparameters h; h.my_beta=beta[b]; h.my_g_thr=1e-6; h.my_maxstep=20000; + CG_Optimizer o; + NRVec x=rosenbrockstart(10); + bool conv; + double f=run(o,&rb,x,h,NULL,&conv); + char s[64]; sprintf(s,"CG variant %d converges on Rosenbrock n=10",b); + check(s,conv && f<1e-10,f); + } +{ +GDIIS_hyperparameters h; h.my_g_thr=1e-6; h.my_maxstep=5000; +GDIIS_Optimizer o; +NRVec x=rosenbrockstart(10); +bool conv; +double f=run(o,&rb,x,h,NULL,&conv); +check("GDIIS converges on Rosenbrock n=10",conv && f<1e-10,f); +//classical GDIIS acceleration with a fixed Hessian (updates switched off) on a well scaled quadratic +class scaledquad : public OptProblem {public: double funcgrad(const NRVec &x, NRVec *g) const {int n=x.size(); double r=0; if(g) {g->copyonwrite(); *g=0.;} + for(int i=0; i a; +x=startpoint(30,1); +f=run(o,&sq,x,h,&a,&conv); +check("GDIIS with a fixed Hessian accelerates (< 100 evaluations on a quadratic n=30)",conv && a.my_maxfunc<100,a.my_maxfunc); +} +{ +//numerical gradients give the same BFGS run as analytic ones +testproblem_vqe vqe; +ParameterShift ps(&vqe); +FiniteDifference fd(&vqe); +BFGS_hyperparameters h; h.my_g_thr=1e-6; +BFGS_Optimizer o; +NRVec x1=startpoint(8,9), x2(x1), x3(x1); +hyperparameters a1, a2, a3; +double f1=run(o,&vqe,x1,h,&a1), f2=run(o,&ps,x2,h,&a2), f3=run(o,&fd,x3,h,&a3); +check("BFGS with ParameterShift and FiniteDifference gradients as with analytic ones",a1.my_maxstep==a2.my_maxstep && a1.my_maxstep==a3.my_maxstep && abs(f1-f2)<1e-12 && abs(f1-f3)<1e-10,f2-f1); +} +} + + +static void test_derivative_free() +{ +cout <<"--- derivative-free methods\n"; +rosenbrock rb; +testproblem tp; +{ +NelderMead_hyperparameters h; h.my_maxstep=100000; +NelderMead_Optimizer o; +NRVec x=rosenbrockstart(2); +bool conv; +double f=run(o,&rb,x,h,NULL,&conv); +check("Nelder-Mead converges on Rosenbrock n=2",conv && f<1e-10,f); +h.my_adaptive=false; +x=rosenbrockstart(2); +f=run(o,&rb,x,h,NULL,&conv); +check("Nelder-Mead with standard coefficients converges on Rosenbrock n=2",conv && f<1e-10,f); +} +{ +Powell_hyperparameters h; +Powell_Optimizer o; +NRVec x=startpoint(10,1); +bool conv; +double f=run(o,&tp,x,h,NULL,&conv); +check("Powell converges on testproblem n=10",conv && f<1e-10,f); +coupledquadratic cq; +x=startpoint(10,1); +f=run(o,&cq,x,h,NULL,&conv); +check("Powell converges on a coupled quadratic n=10",conv && f<1e-10,f); +} +{ +testproblem_vqe vqe; +Rotosolve_hyperparameters h; +Rotosolve_Optimizer o; +NRVec x=startpoint(12,9); +bool conv; +double f=run(o,&vqe,x,h,NULL,&conv); +BFGS_hyperparameters hb; +BFGS_Optimizer ob; +NRVec xb=startpoint(12,9); +double fb=run(ob,&vqe,xb,hb); +check("Rotosolve reaches the minimum of testproblem_vqe found by BFGS",conv && abs(f-fb)<1e-7,f-fb); +trigdegree3 t3; +NRVec y(2); y[0]=0.7; y[1]=0.3; +h.my_maxstep=1; +run(o,&t3,y,h); +double best=1e9; +for(int j=0; j<1000000; ++j) best=min(best,trigdegree3::tr(5.*j/1000000)); +check("Rotosolve degree-3 reconstruction finds the exact minimum along the parameter",abs(trigdegree3::tr(y[0])-best)<1e-10,trigdegree3::tr(y[0])-best); +} +//COBYLA, both radius modes +for(int fl=0; fl<2; ++fl) + { + Cobyla_hyperparameters h; h.my_rhoend=1e-8; h.my_flexible_radius=fl; + Cobyla_Optimizer o; + const char *mode= fl ? "flexible radius" : "original"; + struct {const char *name; const OptProblem *p; NRVec x0; double fmin;} cases[]={ + {"quadratic with a linear constraint",new quadratic_linearconstraint,NRVec(0.,2),2}, + {"x0 x1 in the unit disk",new product_disk,NRVec(1.,2),-0.5}, + {"Rosen-Suzuki",new rosensuzuki,NRVec(0.,4),-44}, + {"Rosen-Suzuki from an infeasible point",new rosensuzuki,NRVec(3.,4),-44}, + {"equality constraint x0^2+x1^2=1",new sum_circle,NRVec(0.5,2),-sqrt(2.)}}; + for(auto &c: cases) + { + NRVec x(c.x0); + bool conv; + double f=run(o,c.p,x,h,NULL,&conv); + const OptConstraints *con=dynamic_cast *>(c.p); + NRVec cv(con->nconstraints()); + con->constraints(x,cv); + string s=string("COBYLA (")+mode+") solves "+c.name; + check(s.c_str(),conv && abs(f-c.fmin)<1e-7 && cv.min()> -1e-10,f-c.fmin); + delete c.p; + } + testproblem_vqe_constrained vc; + NRVec x=startpoint(8,9); + bool conv; + double f=run(o,&vc,x,h,NULL,&conv); + NRVec cv(1); + vc.constraints(x,cv); + string s=string("COBYLA (")+mode+") on testproblem_vqe_constrained n=8"; + check(s.c_str(),conv && abs(f+3.334983606)<1e-7 && cv[0]> -1e-10,f); + } +} + + +static void test_firstorder() +{ +cout <<"--- first-order and stochastic methods\n"; +rosenbrock rb; +NRVec x0=rosenbrockstart(5); +const int steps=300; +//Adam, AdamW, AMSGrad, SGD, Nesterov vs the reference implementation +for(int type=0; type<5; ++type) + { + NRVec x(x0); + double lr= type<3 ? 0.01 : 1e-4; + int rtype= type==2 ? 1 : type==3 ? 2 : type==4 ? 3 : 0; + double wd= type==1 ? 0.1 : 0; + if(type<3) + { + Adam_hyperparameters h; h.my_maxstep=steps; h.my_g_thr=0; h.my_lr=lr; h.my_weight_decay=wd; h.my_amsgrad= type==2; + Adam_Optimizer o; + run(o,&rb,x,h); + } + else + { + SGD_hyperparameters h; h.my_maxstep=steps; h.my_g_thr=0; h.my_lr=lr; h.my_nesterov= type==4; + SGD_Optimizer o; + run(o,&rb,x,h); + } + NRVec xr=reference_firstorder(rb,x0,steps,lr,rtype,wd); + xr-=x; + const char *names[]={"Adam","AdamW","AMSGrad","SGD with momentum","SGD-Nesterov"}; + string s=string(names[type])+" agrees with the reference update formulas"; + check(s.c_str(),xr.norm()<1e-10,xr.norm()); + } +check("x0 unchanged by the in-place updates of copies (copy-on-write)",x0[0]== -1.2 && x0[1]==1); +{ +noisyquadratic nq(0.1); +Adam_hyperparameters h; h.my_maxstep=3000; h.my_g_thr=0; h.my_lr=0.01; h.my_lr_decay=0.01; +Adam_Optimizer o; +NRVec x=startpoint(20,5); +run(o,&nq,x,h); +check("Adam with decaying learning rate reaches the noise floor of a noisy quadratic",noisyquadratic::exact(x)<1e-3,noisyquadratic::exact(x)); +} +{ +//SPSA +noisyquadratic nq(0.1); +SPSA_hyperparameters h; h.my_maxfunc=4000; +SPSA_Optimizer o; +NRVec x=startpoint(10,5); +run(o,&nq,x,h); +check("SPSA on a noisy quadratic n=10",noisyquadratic::exact(x)<0.05,noisyquadratic::exact(x)); +coupledquadratic cq; +double sum=0; +NRVec xs=startpoint(6,3); +for(int r=0; r<100; ++r) + { + SPSA_hyperparameters hc; hc.my_maxstep=1; hc.my_target_step=0.1; + NRVec y(xs); + run(o,&cq,y,hc); + NRVec d=y-xs; + sum+=d.norm(); + } +check("SPSA calibration gives the requested first step length (within 20%)",abs(sum/100-0.1)<0.02,sum/100); +} +{ +//QNG: one step by hand and convergence +testproblem_vqe vqe; +NRVec x=startpoint(8,4,2.), x0q(x); +x0q.copyonwrite(); +QNG_hyperparameters h; h.my_maxstep=1; h.my_g_thr=0; h.my_lr=0.1; h.my_metric_reg=1e-3; +QNG_Optimizer o; +run(o,&vqe,x,h); +NRVec g(8); +vqe.funcgrad(x0q,&g); +NRSMat G(8); +vqe.metric(x0q,G); +G+=1e-3; +linear_solve(G,g); +NRVec xr(x0q); +xr.axpy(-0.1,g); +xr-=x; +check("QNG step equals x - lr (G + lambda I)^-1 g",xr.norm()<1e-14,xr.norm()); +QNG_hyperparameters hc; hc.my_maxstep=10000; +NRVec y=startpoint(12,9); +bool conv; +double f=run(o,&vqe,y,hc,NULL,&conv); +check("QNG converges on testproblem_vqe",conv && f< -5.4,f); +} +} + + +static void test_maxfunc_and_parameters() +{ +cout <<"--- evaluation budget and hyperparameters\n"; +testproblem tp; +testproblem_vqe_constrained vc; +NRVec x=startpoint(6,2), xv=startpoint(8,9); +hyperparameters h; h.my_f_thr=1e-10; h.my_x_thr=1e-6; h.my_maxstep=100000; +BFGS_Optimizer bfgs; LBFGS_Optimizer lbfgs; CG_Optimizer cg; GDIIS_Optimizer gdiis; +NelderMead_Optimizer nm; Powell_Optimizer powell; Adam_Optimizer adam; Rotosolve_Optimizer roto; Cobyla_Optimizer cobyla; +struct {const char *name; Optimizer *o; const OptProblem *p; NRVec *x; int minfunc;} cases[]={ + {"BFGS",&bfgs,&tp,&x,1},{"L-BFGS",&lbfgs,&tp,&x,1},{"CG",&cg,&tp,&x,1},{"GDIIS",&gdiis,&tp,&x,1},{"Nelder-Mead",&nm,&tp,&x,7}, + {"Powell",&powell,&tp,&x,1},{"Adam",&adam,&tp,&x,1},{"Rotosolve",&roto,&vc,&xv,1},{"COBYLA",&cobyla,&vc,&xv,9}}; +for(auto &c: cases) + { + int bad=maxfunc_violations(*c.o,c.p,*c.x,h,c.minfunc,120); + string s=string(c.name)+" respects maxfunc and reports the evaluations"; + check(s.c_str(),bad==0,bad); + } +//parameters given at different levels of the hierarchy, achieved filled at its level +linesearch_hyperparameters lp; lp.my_c2=0.123; lp.my_g_thr=1e-4; +BFGS_hyperparameters a; +NRVec y=startpoint(5,1); +run(bfgs,&tp,y,lp,&a); +check("parameters taken from the caller's level, defaults kept for the rest",a.my_c2==0.123 && a.my_sy_thr==1e-12 && a.my_g_thr<1e-4); +} + + +static void test_diis() +{ +cout <<"--- DIIS\n"; +//DIIS acceleration of a linear fixed-point iteration x = Ax + b +struct fixedpoint {static NRVec run(DIIS,NRVec > &d, const int seed) + { + const int n=5; + srand(seed); + NRMat A(n,n); A.randomize(0.3); + NRVec b(n); b.randomize(1.); + NRVec x(0.,n); + for(int it=0; it<6; ++it) {NRVec xn=A*x; xn+=b; NRVec err=xn-x; d.extrapolate(xn,err); x=xn;} + return x; + }}; +DIIS,NRVec > din(4,true), dout(4,false), dfresh(4,true); +NRVec a=fixedpoint::run(din,1), b=fixedpoint::run(dout,1); +a-=b; +check("DIIS in-core and out-of-core (AuxStorage) give the same extrapolation",a.norm()==0,a.norm()); +din.reset(); +NRVec c=fixedpoint::run(din,2), f=fixedpoint::run(dfresh,2); +c-=f; +check("DIIS after reset() as a fresh object",c.norm()==0,c.norm()); +} + + +//a member function applied to an object sharing its data with a copy must not change the copy +template static bool unchanged(const X &orig, const X &keep) {X d=orig-keep; return d.norm()==0;} + +static void test_copyonwrite() +{ +cout <<"--- copy-on-write in the functions of the vector and matrix classes\n"; +#define FRESHMAT NRMat orig(4,5); orig.randomize(1.); NRMat keep(orig); keep.copyonwrite(); NRMat a(orig); +#define FRESHVEC NRVec orig(6); orig.randomize(1.); NRVec keep(orig); keep.copyonwrite(); NRVec a(orig); +#define FRESHSMAT NRSMat orig(5); orig.randomize(1.); NRSMat keep(orig); keep.copyonwrite(); NRSMat a(orig); +NRVec r(5); r.randomize(1.); +NRVec col(4); col.randomize(1.); +{FRESHMAT a.rowset(r,1); check("NRMat::rowset",unchanged(orig,keep));} +{FRESHMAT a.columnset(col,2); check("NRMat::columnset",unchanged(orig,keep));} +{FRESHMAT NRMat s(2,2); s=1.; a.storesubmatrix(1,1,s); check("NRMat::storesubmatrix(row,col)",unchanged(orig,keep));} +{FRESHMAT NRMat s(2,2); s=1.; NRVec ri(2), ci(2); ri[0]=0; ri[1]=3; ci[0]=1; ci[1]=4; a.storesubmatrix(ri,ci,s); check("NRMat::storesubmatrix(rows,cols)",unchanged(orig,keep));} +{FRESHMAT a.randomize(1.); check("NRMat::randomize",unchanged(orig,keep));} +{NRMat orig(4,4); orig.randomize(1.); NRMat keep(orig); keep.copyonwrite(); NRMat a(orig); NRPerm p(4); p.identity(); a.axpy(2.,p,true); check("NRMat::axpy(permutation)",unchanged(orig,keep));} +{NRMat orig(4,4); orig.randomize(1.); NRMat keep(orig); keep.copyonwrite(); NRMat a(orig); a.orthonormalize(true,NULL); check("NRMat::orthonormalize",unchanged(orig,keep));} +{ +NRMat m(4,4); m.randomize(1.); +NRMat t(m); t.copyonwrite(); t.transposeme(3); +bool ok=true; +for(int i=0; i<4; ++i) for(int j=0; j<4; ++j) if(t(i,j)!=((i<3 && j<3) ? m(j,i) : m(i,j))) ok=false; +check("NRMat::transposeme(n) transposes the leading n x n block",ok); +} +{ +FILE *f=tmpfile(); +fprintf(f,"2 3\n1 2 3\n4 5 6\n"); +rewind(f); +NRMat m; +m.fscanf(f,"%lf"); +fclose(f); +check("NRMat::fscanf reads a non-square matrix",m.nrows()==2 && m.ncols()==3 && m(0,2)==3 && m(1,0)==4 && m(1,2)==6); +} +{FRESHVEC NRVec s(2); s=1.; a.storesubvector(1,s); check("NRVec::storesubvector(from)",unchanged(orig,keep));} +{FRESHVEC NRVec s(2); s=1.; NRVec sel(2); sel[0]=0; sel[1]=4; a.storesubvector(sel,s); check("NRVec::storesubvector(selection)",unchanged(orig,keep));} +{FRESHVEC a.randomize(1.); check("NRVec::randomize",unchanged(orig,keep));} +{FRESHVEC NRVec w(6); w.randomize(1.); a.call_on_me([](const double &p, const double &q){return p>q ? p : q;},w); check("NRVec::call_on_me binary",unchanged(orig,keep));} +{FRESHSMAT a.randomize(1.); check("NRSMat::randomize",unchanged(orig,keep));} +{FRESHSMAT NRSMat s(2); s=1.; a.storesubmatrix(1,s); check("NRSMat::storesubmatrix(from)",unchanged(orig,keep));} +{ +FRESHSMAT +NRSMat s(3); s.randomize(1.); +NRVec sel(3); sel[0]=0; sel[1]=2; sel[2]=4; +a.storesubmatrix(sel,s); +bool ok=true; +for(int i=0; i<3; ++i) for(int j=0; j<3; ++j) if(a(sel[i],sel[j])!=s(i,j)) ok=false; +check("NRSMat::storesubmatrix(selection) copy-on-write and size check",unchanged(orig,keep) && ok); +} +#undef FRESHMAT +#undef FRESHVEC +#undef FRESHSMAT +} + + +int main(int argc, char **argv) +{ +test_problems_and_capabilities(); +test_linemin(); +test_gradient_methods(); +test_derivative_free(); +test_firstorder(); +test_maxfunc_and_parameters(); +test_diis(); +test_copyonwrite(); +cout <<"=== "< or + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +//test problems for the optimization methods, used by t.cc and testoptimize.cc + +#ifndef _TESTOPTIMIZE_H_ +#define _TESTOPTIMIZE_H_ + +#include "optimize.h" + +namespace LA { + +//test case +template +class testproblem : public OptProblem { +public: +virtual T funcgrad(const NRVec &x, NRVec *grad) const; +virtual ~testproblem() {}; +}; + +//test case with all capabilities: Ising chain with longitudinal and transverse fields in the product state of qubits RY(b_i) RX(a_i)|0>, +//x = (a_0,b_0,a_1,b_1,...), E = sum_i (h_i z_i + k_i x_i) + sum_i J_i z_i z_{i+1}, Bloch vector components z_i = cos a_i cos b_i, x_i = cos a_i sin b_i; +//trigonometric of degree 1 with period 2 pi in all parameters, Fubini-Study metric 1/4 diag(1, cos^2 a_i) for each qubit +template +class testproblem_vqe : public OptProblem, public OptHessian, public OptMetric, public OptTrigonometric { +public: +static T hfield(const int i) {return 0.7+0.1*(i%3);} +static T kfield(const int i) {return 0.4;} +static T J(const int i) {return 1;} +virtual T funcgrad(const NRVec &x, NRVec *grad) const; +virtual void hessian(const NRVec &x, NRSMat &h) const; +virtual void metric(const NRVec &x, NRSMat &g) const; +virtual int degree(const int i) const {return 1;} +virtual ~testproblem_vqe() {}; +}; + +//constrained test case: testproblem_vqe with nonnegative total magnetization sum_i z_i >= 0 +template +class testproblem_vqe_constrained : public testproblem_vqe, public OptConstraints { +public: +virtual int nconstraints() const {return 1;} +virtual void constraints(const NRVec &x, NRVec &c) const; +virtual ~testproblem_vqe_constrained() {}; +}; + + +template +T testproblem::funcgrad(const NRVec &x, NRVec *grad) const +{ +int n=x.size(); +T r= n*(n-1)/2; +for(int i=0; isize()!=x.size()) laerror("grad size mismatch"); + grad->copyonwrite(); //I forgot this initially + for(int k=0; k +void testproblem_vqe_constrained::constraints(const NRVec &x, NRVec &c) const +{ +c.copyonwrite(); +c[0]=0; +for(int i=0; i and their derivatives +template +class blochqubit { +public: +T z, x, za, zb, xa, xb; +blochqubit(const T a, const T b) + { + T ca=std::cos(a), sa=std::sin(a), cb=std::cos(b), sb=std::sin(b); + z=ca*cb; x=ca*sb; za= -sa*cb; zb= -ca*sb; xa= -sa*sb; xb=ca*cb; + } +}; + +template +T testproblem_vqe::funcgrad(const NRVec &x, NRVec *grad) const +{ +const int nq=x.size()/2; +if(2*nq!=x.size()) laerror("testproblem_vqe needs an even number of parameters"); +T r=0; +for(int i=0; i q(x[2*i],x[2*i+1]); + r+= hfield(i)*q.z + kfield(i)*q.x; + if(isize()!=x.size()) laerror("grad size mismatch"); + grad->copyonwrite(); + for(int i=0; i q(x[2*i],x[2*i+1]); + T eff=hfield(i); //coefficient of z_i including the couplings + if(i>0) eff+= J(i-1)*std::cos(x[2*i-2])*std::cos(x[2*i-1]); + if(i +void testproblem_vqe::hessian(const NRVec &x, NRSMat &h) const +{ +const int nq=x.size()/2; +h=(T)0; +for(int i=0; i q(x[2*i],x[2*i+1]); + T eff=hfield(i); + if(i>0) eff+= J(i-1)*std::cos(x[2*i-2])*std::cos(x[2*i-1]); + if(i q1(x[2*i+2],x[2*i+3]); + h(2*i,2*i+2)= J(i)*q.za*q1.za; + h(2*i,2*i+3)= J(i)*q.za*q1.zb; + h(2*i+1,2*i+2)= J(i)*q.zb*q1.za; + h(2*i+1,2*i+3)= J(i)*q.zb*q1.zb; + } + } +} + +template +void testproblem_vqe::metric(const NRVec &x, NRSMat &g) const +{ +const int nq=x.size()/2; +g=(T)0; +for(int i=0; i