From 3ef603c1840d3965096169c2f7a07a72bc208bbb Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Mon, 26 Oct 2026 16:31:36 +0100 Subject: [PATCH] optimize: freezing parameters in a problem using a bitvector --- optimize.cc | 60 +++++++++++++++++++++++++++++++++++ optimize.h | 33 ++++++++++++++++++-- t.cc | 22 +++++++++++++ testoptimize.cc | 83 +++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 196 insertions(+), 2 deletions(-) diff --git a/optimize.cc b/optimize.cc index f17ce6a..8c300e9 100644 --- a/optimize.cc +++ b/optimize.cc @@ -24,6 +24,65 @@ namespace LA { +template +void FrozenOptProblem::setfrozen(const bitvector &frozen) +{ +const int n=xfull.size(); +if((int)frozen.size()!=n) laerror("size mismatch in FrozenOptProblem::setfrozen"); +NRVec index(n-(int)frozen.population()); +for(int i=0, k=0; i +T FrozenOptProblem::funcgrad(const NRVec &x, NRVec *grad) const +{ +NRVec full=expand(x); +if(!grad) return this->p->funcgrad(full,NULL); +NRVec gfull(full.size()); +T f=this->p->funcgrad(full,&gfull); +grad->subvector_from(freeindex,gfull); +return f; +} + +//the capabilities are used only if the wrapped problem has them, as ensured by optcapability() +template +void FrozenOptProblem::hessian(const NRVec &x, NRSMat &h) const +{ +NRSMat hfull(xfull.size()); +optcapability >(this->p)->hessian(expand(x),hfull); +h=hfull.submatrix(freeindex); +} + +template +NRVec FrozenOptProblem::hessian_times(const NRVec &x, const NRVec &v) const +{ +NRVec vfull((T)0,xfull.size()); +v.storesubvector_to(freeindex,vfull); +return optcapability >(this->p)->hessian_times(expand(x),vfull).subvector(freeindex); +} + +template +void FrozenOptProblem::metric(const NRVec &x, NRSMat &g) const +{ +NRSMat gfull(xfull.size()); +optcapability >(this->p)->metric(expand(x),gfull); +g=gfull.submatrix(freeindex); +} + +template +int FrozenOptProblem::degree(const int i) const {return optcapability >(this->p)->degree(freeindex[i]);} + +template +T FrozenOptProblem::period(const int i) const {return optcapability >(this->p)->period(freeindex[i]);} + +template +int FrozenOptProblem::nconstraints() const {return optcapability >(this->p)->nconstraints();} + +template +void FrozenOptProblem::constraints(const NRVec &x, NRVec &c) const {optcapability >(this->p)->constraints(expand(x),c);} + + //central difference (f(x+h e_i) - f(x-h e_i))/(2h); xs is a working copy of x, restored on return template static T central_difference(const OptProblem *p, NRVec &xs, const NRVec &x, const int i, const T h) @@ -1553,6 +1612,7 @@ template double linesearch_wolfe(optevaluator &e, NRVec &x, doub template class optstate; template class ParameterShift; template class FiniteDifference; +template class FrozenOptProblem; template class BFGS_Optimizer; template class LBFGS_Optimizer; template class GDIIS_Optimizer; diff --git a/optimize.h b/optimize.h index ef85e8c..4f8e93d 100644 --- a/optimize.h +++ b/optimize.h @@ -28,6 +28,7 @@ #include "vec.h" #include "mat.h" #include "smat.h" +#include "bitvector.h" namespace LA { @@ -90,19 +91,47 @@ class OptProblemAdapter : public OptProblem { public: const OptProblem *p; //wrapped problem OptProblemAdapter(const OptProblem *pp) : p(pp) {}; +virtual bool transforms_capabilities() const {return false;} //the adapter provides the capabilities of the wrapped problem in a transformed form (e.g. other dimension), not the original ones }; //capability C (OptHessian, OptMetric, ...) of the problem, searched also through adapters; NULL if not available +//a transforming adapter provides C only if the wrapped problem does, and never passes on the untransformed one template const C *optcapability(const OptProblem *o) { const C *c=dynamic_cast(o); -if(c) return c; const OptProblemAdapter *a=dynamic_cast *>(o); -return a ? optcapability(a->p) : NULL; +if(!a) return c; +if(a->transforms_capabilities()) return c && optcapability(a->p) ? c : NULL; +return c ? c : optcapability(a->p); } +//optimization in the free parameters only, the frozen ones keep their values from the full parameter vector +//the Hessian, metric, trigonometric structure and constraints of the wrapped problem are provided for the free parameters +//usage: FrozenOptProblem fp(&problem,x,frozen); NRVec xr=fp.reduce(x); optimizer.minimize(&fp,xr); x=fp.expand(xr); +template +class FrozenOptProblem : public OptProblemAdapter, public OptHessian, public OptMetric, public OptTrigonometric, public OptConstraints { +public: +NRVec xfull; //full parameter vector, the values of the frozen parameters are taken from it +NRVec freeindex; //indices of the free parameters in the full vector +FrozenOptProblem(const OptProblem *pp, const NRVec &x, const bitvector &frozen) : OptProblemAdapter(pp), xfull(x) {setfrozen(frozen);}; +void setfrozen(const bitvector &frozen); //set of frozen parameters, bit i set = x_i frozen +void setfull(const NRVec &x) {xfull=x;} //new values of the frozen parameters +NRVec reduce(const NRVec &full) const {return full.subvector(freeindex);} //free parameters of a full vector +NRVec expand(const NRVec &reduced) const {NRVec full(xfull); reduced.storesubvector_to(freeindex,full); return full;} //full vector with the frozen values +virtual bool transforms_capabilities() const {return true;} +virtual T funcgrad(const NRVec &x, NRVec *grad) const; +virtual void hessian(const NRVec &x, NRSMat &h) const; +virtual NRVec hessian_times(const NRVec &x, const NRVec &v) const; +virtual void metric(const NRVec &x, NRSMat &g) const; +virtual int degree(const int i) const; +virtual T period(const int i) const; +virtual int nconstraints() const; +virtual void constraints(const NRVec &x, NRVec &c) const; +}; + + //gradient by the general parameter-shift rule (Wierichs et al., Quantum 6, 677 (2022); degree 1: Mitarai et al., Phys. Rev. A 98, 032309 (2018), //Schuld et al., Phys. Rev. A 99, 032331 (2019)), exact from 2R_i evaluations per parameter; parameters of degree 0 by central differences with step h //the wrapped problem has to provide OptTrigonometric diff --git a/t.cc b/t.cc index f9311aa..108619a 100644 --- a/t.cc +++ b/t.cc @@ -4882,6 +4882,28 @@ double r=problem.funcgrad(x,NULL); cout < x(n); +x.randomize(1.); +testproblem problem; +bitvector frozen(n); +for(int i=0; i fproblem(&problem,x,frozen); +NRVec xr=fproblem.reduce(x); +BFGS_hyperparameters param; +param.my_g_thr=1e-6; +cout << "initial " < o; +optresult res=o.minimize(&fproblem,xr,¶m); +NRVec xopt=fproblem.expand(xr); +cout << "optimized " < vqe; +const int n=8; +NRVec x=startpoint(n,4,2.); +bitvector frozen(n); +frozen.set(1); frozen.set(4); frozen.set(5); +FrozenOptProblem fp(&vqe,x,frozen); +NRVec xr=fp.reduce(x); +NRVec xe=fp.expand(xr); +xe-=x; +check("reduce and expand",xr.size()==5 && xe.norm()==0); +NRVec y=startpoint(5,7); //other values of the free parameters +NRVec yfull=fp.expand(y); +NRVec gfull(n), g(5); +double ffull=vqe.funcgrad(yfull,&gfull), f=fp.funcgrad(y,&g); +NRVec gd=gfull.subvector(fp.freeindex)-g; +check("reduced function value and gradient",f==ffull && gd.norm()==0); +NRSMat hfull(n), h(5), mfull(n), m(5); +vqe.hessian(yfull,hfull); fp.hessian(y,h); +vqe.metric(yfull,mfull); fp.metric(y,m); +NRSMat hd=hfull.submatrix(fp.freeindex)-h, md=mfull.submatrix(fp.freeindex)-m; +check("reduced Hessian and metric",hd.norm()==0 && md.norm()==0); +NRVec v=startpoint(5,8), hv=fp.hessian_times(y,v), hv2=h*v; +hv-=hv2; +check("reduced Hessian-vector product",hv.norm()<1e-14,hv.norm()); +NRVec cfull(1), c(1); +vqe.constraints(yfull,cfull); fp.constraints(y,c); +check("constraints and trigonometric degree of the reduced problem",fp.nconstraints()==1 && c[0]==cfull[0] && fp.degree(2)==1 && fp.period(2)==2*M_PI); +testproblem tp; +FrozenOptProblem ftp(&tp,startpoint(6,1),bitvector(6)); +ParameterShift ps(&fp); +check("capabilities of the frozen problem only if the wrapped problem has them",optcapability >(&fp)==&fp && optcapability >(&ftp)==NULL && optcapability >(&ftp)==NULL && optcapability >(&ps)==&fp); +//optimizations keep the frozen parameters +{ +bitvector fr(10); +for(int i=0; i<10; i+=3) fr.set(i); +NRVec x0=startpoint(10,1); +FrozenOptProblem fq(&tp,x0,fr); +NRVec z=fq.reduce(x0); +BFGS_hyperparameters hb; hb.my_g_thr=1e-8; +BFGS_Optimizer o; +bool conv; +run(o,&fq,z,hb,NULL,&conv); +NRVec zf=fq.expand(z), gz(z.size()); +fq.funcgrad(z,&gz); +bool kept=true; +const NRVec &x0c=x0; //x0 is shared with fq.xfull, read-only access +for(int i=0; i<10; i+=3) if(zf[i]!=x0c[i]) kept=false; +check("BFGS with frozen parameters: frozen ones kept, reduced gradient zero",conv && kept && gz.norm()<1e-8,gz.norm()); +} +{ +NRVec z=fp.reduce(x); +QNG_hyperparameters hq; hq.my_maxstep=10000; +QNG_Optimizer o; +bool conv; +testproblem_vqe vu; +FrozenOptProblem fu(&vu,x,frozen); +run(o,&fu,z,hq,NULL,&conv); +NRVec gz(5); +fu.funcgrad(z,&gz); +check("QNG uses the reduced metric",conv && gz.norm()<1e-5,gz.norm()); +z=fp.reduce(x); +Rotosolve_hyperparameters hr; +Rotosolve_Optimizer orot; +run(orot,&fu,z,hr,NULL,&conv); +fu.funcgrad(z,&gz); +check("Rotosolve uses the reduced trigonometric structure",conv && gz.norm()<1e-4,gz.norm()); +z=fp.reduce(x); +Cobyla_hyperparameters hc; hc.my_rhoend=1e-8; +Cobyla_Optimizer oc; +run(oc,&fp,z,hc,NULL,&conv); +NRVec cz(1); +fp.constraints(z,cz); +NRVec zf=fp.expand(z); +const NRVec &xc=x; //x is shared with fp.xfull, read-only access +check("COBYLA with frozen parameters satisfies the constraint and keeps them",conv && cz[0]> -1e-10 && zf[1]==xc[1] && zf[4]==xc[4] && zf[5]==xc[5],cz[0]); +} +} + + static void test_maxfunc_and_parameters() { cout <<"--- evaluation budget and hyperparameters\n"; @@ -725,6 +807,7 @@ test_linemin(); test_gradient_methods(); test_derivative_free(); test_firstorder(); +test_frozen(); test_maxfunc_and_parameters(); test_diis(); test_copyonwrite();