From 2e6ce0c13e3e125060caecbec9c2cf256e997169 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Mon, 26 Oct 2026 15:29:48 +0100 Subject: [PATCH] optimize: cobyla vibe-coded --- optimize.cc | 323 +++++++++++++++++++++++++++++++++++++++++++++++++++- optimize.h | 43 +++++++ t.cc | 9 +- 3 files changed, 370 insertions(+), 5 deletions(-) diff --git a/optimize.cc b/optimize.cc index 4717153..aa3f599 100644 --- a/optimize.cc +++ b/optimize.cc @@ -49,6 +49,15 @@ return r; } +template +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) @@ -156,7 +165,6 @@ template void testproblem_vqe::hessian(const NRVec &x, NRSMat &h) const { const int nq=x.size()/2; -h.copyonwrite(); h=(T)0; for(int i=0; i void testproblem_vqe::metric(const NRVec &x, NRSMat &g) const { const int nq=x.size()/2; -g.copyonwrite(); g=(T)0; for(int i=0; i= r and |(z_0,...,z_{nd-1})| <= delta, nd = nz or nz-1, from a feasible z by an active-set method: +//z moves along the projection of -c orthogonal to the normals of the active constraints until a further constraint becomes active, +//the trust region boundary is reached (which finishes the calculation as in COBYLA), or the multipliers show optimality; +//the projection is computed by least squares, which tolerates dependent normals at degenerate points +template +static void lp_trustregion(const NRVec &c, const NRMat &A, const NRVec &r, const int nd, const T delta, NRVec &z) +{ +const int nz=c.size(); +const int m=r.size(); +const T tiny=1e3*std::numeric_limits::epsilon()*c.norm(); +NRVec active(m); +active.clear(); +NRVec cols(nz); +for(int j=0; j &zc=z; //element access for reading only +const int maxit=2*(m+nz)+10; +for(int it=0; it rows(k); + k=0; + for(int i=0; i s = c*(T)-1; //search direction + NRVec lambda; //multipliers of the active constraints + if(k>0) + { + NRMat N=A.submatrix(rows,cols); //normals of the active constraints + NRMat Nt=N.transpose(); + lambda=svd_solve(Nt,c); //least squares c ~ N^T lambda + NRVec nl=lambda*N; + s+=nl; + } + if(s.norm()<=tiny) //stationary on the active set, optimal if the multipliers are nonnegative + { + int drop= -1; + T lmin=0; + for(int j=0; j Az=A*z; + NRVec As=A*s; + T alpha=std::numeric_limits::max(); + int block= -1; + for(int i=0; i0) + { + T ab=(-zs+std::sqrt(std::max((T)0,zs*zs-ss*(zz-delta*delta))))/ss; + if(ab= 0 within |d| <= rho; if they cannot be satisfied +//in the trust region, the largest violation t is minimized first (stage 1) and then kept, c + A d >= -t (stage 2) (Powell 1994) +template +static NRVec cobyla_step(const NRVec &g, const NRMat &A, const NRVec &c, const T rho) +{ +const int n=g.size(); +const int m=c.size(); +NRVec d((T)0,n); +if(m==0) {T gn=g.norm(); if(gn>0) d.axpy(-rho/gn,g); return d;} +T t=std::max((T)0,-c.min()); +if(t>0) //stage 1 in the variables (d,t): min t subject to c_i + a_i.d + t >= 0, t >= 0 + { + NRMat A1(m+1,n+1); + A1.clear(); + A1.storesubmatrix(0,0,A); + for(int i=0; i<=m; ++i) A1(i,n)=1; + NRVec r1(m+1); //r1_i = -c_i, r1_m = 0 + r1.clear(); + for(int i=0; i c1((T)0,n+1); + c1[n]=1; + NRVec z((T)0,n+1); + z[n]=t; + lp_trustregion(c1,A1,r1,n,rho,z); + t=z[n]; + for(int j=0; j= -t +NRVec r2 = c*(T)-1; +r2-=t; +lp_trustregion(g,A,r2,n,rho,d); +return d; +} + + +template +static T cobyla_violation(const NRVec &c) {return c.size()>0 ? std::max((T)0,-c.min()) : 0;} + +//f and the constraint values at x +template +static T cobyla_eval(optevaluator &e, const OptConstraints *con, const NRVec &x, NRVec &c) +{ +T f=e(x); +if(con) con->constraints(x,c); +return f; +} + + +template +optresult Cobyla_Optimizer::minimize(const OptProblem *o, NRVec &x, const hyperparameters *par, hyperparameters *achieved) +{ +const char *name="COBYLA"; +Cobyla_hyperparameters hp; +if(par) hp.setfrom(*par); +if(hp.my_saddle_order) laerror("COBYLA cannot search for saddle points"); +const T acc_sigma=0.25, acc_eta=2.1, geo_gamma=0.5; //simplex acceptability and geometry step (Powell 1994) +const int verbose=hp.my_verbose; +const int n=x.size(); +const OptConstraints *con=optcapability >(o); +const int m= con ? con->nconstraints() : 0; +optevaluator e(o,hp.my_maxfunc); +optstate st(name,hp,e); + +//initial simplex x, x + rhobeg e_j with function and constraint values +T rho=hp.my_rhobeg; //resolution, lower bound of the trust region radius +T delta=rho; //trust region radius, always equal to rho in Powell's original COBYLA +const bool flexible=hp.my_flexible_radius; +T mu=0; //penalty parameter of the merit function f + mu*(largest constraint violation) +NRVec > v(n+1), cv(n+1); +NRVec fv(n+1); +for(int j=0; j<=n; ++j) + { + v[j]=x; + if(j>0) {v[j].copyonwrite(); v[j][j-1]+=rho;} + cv[j]=NRVec(m); + fv[j]=cobyla_eval(e,con,v[j],cv[j]); + } + +int best=0; +T fprev=fv[0]; +NRVec other(n); //vertices other than the best one +bool lastsuccess=false; //the last trust region step was successful, the next one follows without checking the simplex +for(;; ++st.iter) + { + best=0; + for(int j=1; j<=n; ++j) if(fv[j]+mu*cobyla_violation(cv[j]) < fv[best]+mu*cobyla_violation(cv[best])) best=j; + if(st.iter>0) {st.step(fprev-fv[best],rho); fprev=fv[best];} + if(st.stop(fv[best])) break; + for(int j=0, k=0; j<=n; ++j) if(j!=best) other[k++]=j; + + //displacements d_k of the other vertices from the best one and the dual basis B = D^-T, b_k.d_l = delta_kl + NRMat D(n,n); + for(int k=0; k dk=v[other[k]]-v[best]; D.rowset(dk,k);} + NRMat B=D.inverse().transpose(); + //linear models: gradient g of f, gradients of the constraints in the rows of A, values c0 at the best vertex + NRVec df(n); + for(int k=0; k g=df*B; + const NRVec c0(cv[best]); + NRMat A; + if(m>0) + { + NRMat dC(n,m); + for(int k=0; k dc=cv[other[k]]-c0; dC.rowset(dc,k);} + A=dC.transpose()*B; + } + + //acceptability of the simplex: distances eta_k = |d_k| <= acc_eta delta and sigma_k = 1/|b_k| >= acc_sigma rho of the vertices from the opposite faces + int worst= -1; + T etamax=acc_eta*delta; + for(int k=0; ketamax) {etamax=eta; worst=k;}} + if(worst<0) + { + T sigmin=acc_sigma*rho; + for(int k=0; k=0) //geometry step: the worst vertex moves to the best one + a step along the normal of its opposite face, with the sign preferred by the merit model + { + NRVec dir=B.row(worst); + //step length geo_gamma rho, or with flexible radius max(min(0.1 DIST, 0.5 delta), rho) (NEWUOA eq. (6.16)), DIST the distance of the replaced vertex + T len= flexible ? std::max(std::min((T)0.1*D.row(worst).norm(),(T)0.5*delta),rho) : geo_gamma*rho; + dir*= len/dir.norm(); + T mplus=g*dir, mminus= -mplus; + if(m>0) + { + NRVec ad=A*dir; + NRVec cp(c0), cm(c0); + cp+=ad; cm-=ad; + mplus+=mu*cobyla_violation(cp); + mminus+=mu*cobyla_violation(cm); + } + if(mminus1) std::cout < xn=v[best]+dir; + NRVec cn(m); + T fn=cobyla_eval(e,con,xn,cn); + v[other[worst]]=xn; fv[other[worst]]=fn; cv[other[worst]]=cn; + continue; + } + + //trust region step + lastsuccess=false; + NRVec d=cobyla_step(g,A,c0,delta); + const T dnorm=d.norm(); + const T deltaold=delta; //trust region radius of this step + T ratio= -1; //reduction ratio of the merit function + const bool shortstep= dnorm<0.5*rho; + bool failed=shortstep; + if(shortstep && flexible) delta=std::max((T)0.1*delta,rho); //NEWUOA Box 15 + if(!failed) + { + const T v0=cobyla_violation(c0); + NRVec cd(c0); + if(m>0) {NRVec ad=A*d; cd+=ad;} + T fpred= -(g*d); //predicted reductions of f and of the violation + T vpred=v0-cobyla_violation(cd); + if(vpred>0 && mu < -1.5*fpred/vpred) //penalty has to make the predicted reduction of the merit function positive + { + mu= -2*fpred/vpred; + if(verbose>1) std::cout < xn=v[best]+d; + NRVec cn(m); + T fn=cobyla_eval(e,con,xn,cn); + T phi0=fv[best]+mu*v0; + T phin=fn+mu*cobyla_violation(cn); + ratio=(phi0-phin)/pred; + if(verbose>1) std::cout < class trigpolynomial : public OptProblem { @@ -1351,6 +1668,7 @@ template double linesearch_wolfe(optevaluator &e, NRVec &x, doub 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; @@ -1360,6 +1678,7 @@ template class CG_Optimizer; template class NelderMead_Optimizer; template class Powell_Optimizer; template class Rotosolve_Optimizer; +template class Cobyla_Optimizer; template class SGD_Optimizer; template class QNG_Optimizer; template class Adam_Optimizer; diff --git a/optimize.h b/optimize.h index 15892bf..33fbf0a 100644 --- a/optimize.h +++ b/optimize.h @@ -74,6 +74,16 @@ virtual ~OptTrigonometric() {}; }; +//inequality constraints c_i(x) >= 0, i=0,...,nconstraints()-1; an equality constraint can be given as two inequalities +template +class OptConstraints { +public: +virtual int nconstraints() const =0; +virtual void constraints(const NRVec &x, NRVec &c) const =0; //c is allocated by the caller with nconstraints() elements +virtual ~OptConstraints() {}; +}; + + //problem wrapping another one (e.g. providing gradients); the capabilities of the wrapped problem remain accessible by optcapability() template class OptProblemAdapter : public OptProblem { @@ -139,6 +149,15 @@ 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 @@ -462,6 +481,30 @@ virtual ~Powell_Optimizer() {}; }; +//COBYLA, constrained optimization by linear approximations (Powell, in Advances in Optimization and Numerical Analysis, eds. Gomez & Hennart, +//Kluwer 1994, p. 51; Powell, Acta Numerica 7, 287 (1998)), derivative-free; constraints from OptConstraints, without them unconstrained +//linear models of f and the constraints are interpolated on a simplex of n+1 points, the trust region radius rho decreases from rhobeg to rhoend; +//the iterations finish when rho reaches rhoend and the trust region step does not improve the point, or by the general criterion +//with f_thr and x_thr, which are used for the change of f in an iteration and for rho (both 0 by default) +template +class Cobyla_hyperparameters : public hyperparameters { +public: +T my_rhobeg; //initial trust region radius, about one tenth of the expected change of x +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;}; +virtual void setfrom(const hyperparameters &p) {const Cobyla_hyperparameters *q=dynamic_cast(&p); if(q) *this= *q; else hyperparameters::setfrom(p);} +}; + +template +class Cobyla_Optimizer : public Optimizer { +public: +virtual optresult minimize(const OptProblem *o, NRVec &x, const hyperparameters *par=NULL, hyperparameters *achieved=NULL); +virtual ~Cobyla_Optimizer() {}; +}; + + //Rotosolve: sweeps of exact minimizations along single parameters in which f is trigonometric (OptTrigonometric) //(Ostaszewski et al., Quantum 5, 391 (2021); Nakanishi et al., Phys. Rev. Research 2, 043158 (2020); degree > 1: Wierichs et al., Quantum 6, 677 (2022)); //f along x_i is reconstructed from 2R_i+1 equidistant samples over its period (2R_i new evaluations) and minimized without further evaluations; diff --git a/t.cc b/t.cc index fe62ec2..de77c6a 100644 --- a/t.cc +++ b/t.cc @@ -4893,17 +4893,20 @@ int n; cin >>n; NRVec x(n); x.randomize(1.); -testproblem_vqe problem; //needs an even dimension +//testproblem_vqe problem; //needs an even dimension +testproblem_vqe_constrained problem; //with the constraint sum_i z_i >= 0 (OptConstraints), for COBYLA ParameterShift psproblem(&problem); //gradient by the parameter-shift rule from function values only //QNG_hyperparameters param; //tested -Rotosolve_hyperparameters param; +//Rotosolve_hyperparameters param; //tested +Cobyla_hyperparameters param; param.my_maxstep=2000; param.my_f_thr=1e-6; param.my_x_thr=1e-3; //param.my_g_thr=1e-3; //gradient methods cout << "initial " < o; //tested -Rotosolve_Optimizer o; +//Rotosolve_Optimizer o; //tested +Cobyla_Optimizer o; //gradient_hyperparameters achieved; //gradient methods hyperparameters achieved; //derivative-free methods optresult res=o.minimize(&problem,x,¶m,&achieved);