diff --git a/optimize.cc b/optimize.cc
index cdd3cfc..4717153 100644
--- a/optimize.cc
+++ b/optimize.cc
@@ -18,6 +18,8 @@
 
 #include "optimize.h"
 #include "la_random.h"
+#include "diis.h"
+#include "nonclass.h"
 
 namespace LA {
 
@@ -47,6 +49,151 @@ return r;
 }
 
 
+//central difference (f(x+h e_i) - f(x-h e_i))/(2h); xs is a working copy of x, restored on return
+template<typename T>
+static T central_difference(const OptProblem<T> *p, NRVec<T> &xs, const NRVec<T> &x, const int i, const T h)
+{
+xs[i]=x[i]+h;
+T d=p->funcgrad(xs,NULL);
+xs[i]=x[i]-h;
+d-=p->funcgrad(xs,NULL);
+xs[i]=x[i];
+return d/(2*h);
+}
+
+
+template<typename T>
+ParameterShift<T>::ParameterShift(const OptProblem<T> *pp, const T hh) : OptProblemAdapter<T>(pp), trig(optcapability<OptTrigonometric<T> >(pp)), h(hh)
+{
+if(!trig) laerror("ParameterShift needs a problem providing OptTrigonometric");
+}
+
+//derivative of a trigonometric polynomial of degree R with frequencies 1..R (in units of w):
+//f'(0) = sum_{mu=1}^{2R} (-1)^{mu-1} f(s_mu)/(4R sin^2(s_mu/2)), s_mu = (2mu-1)pi/(2R) (Wierichs et al., Quantum 6, 677 (2022))
+template<typename T>
+T ParameterShift<T>::funcgrad(const NRVec<T> &x, NRVec<T> *grad) const
+{
+T f=this->p->funcgrad(x,NULL);
+if(!grad) return f;
+grad->copyonwrite();
+NRVec<T> xs(x);
+xs.copyonwrite(); //working copy with shifted elements
+for(int i=0; i<x.size(); ++i)
+	{
+	const int R=trig->degree(i);
+	if(R==0) {(*grad)[i]=central_difference(this->p,xs,x,i,h); continue;}
+	const T w=2*M_PI/trig->period(i);
+	T d=0;
+	for(int mu=1; mu<=2*R; ++mu)
+		{
+		const T s=(2*mu-1)*M_PI/(2*R);
+		const T sn=std::sin(s/2);
+		xs[i]=x[i]+s/w;
+		d+= (mu%2 ? 1 : -1)*this->p->funcgrad(xs,NULL)/(4*R*sn*sn);
+		}
+	xs[i]=x[i];
+	(*grad)[i]=w*d;
+	}
+return f;
+}
+
+
+template<typename T>
+T FiniteDifference<T>::funcgrad(const NRVec<T> &x, NRVec<T> *grad) const
+{
+T f=this->p->funcgrad(x,NULL);
+if(!grad) return f;
+grad->copyonwrite();
+NRVec<T> xs(x);
+xs.copyonwrite(); //working copy with shifted elements
+for(int i=0; i<x.size(); ++i) (*grad)[i]=central_difference(this->p,xs,x,i,h);
+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<typename T>
+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<typename T>
+T testproblem_vqe<T>::funcgrad(const NRVec<T> &x, NRVec<T> *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<nq; ++i)
+	{
+	blochqubit<T> q(x[2*i],x[2*i+1]);
+	r+= hfield(i)*q.z + kfield(i)*q.x;
+	if(i<nq-1) r+= J(i)*q.z*std::cos(x[2*i+2])*std::cos(x[2*i+3]);
+	}
+if(grad)
+	{
+	if(grad->size()!=x.size()) laerror("grad size mismatch");
+	grad->copyonwrite();
+	for(int i=0; i<nq; ++i)
+		{
+		blochqubit<T> 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<nq-1) eff+= J(i)*std::cos(x[2*i+2])*std::cos(x[2*i+3]);
+		(*grad)[2*i]= eff*q.za + kfield(i)*q.xa;
+		(*grad)[2*i+1]= eff*q.zb + kfield(i)*q.xb;
+		}
+	}
+return r;
+}
+
+template<typename T>
+void testproblem_vqe<T>::hessian(const NRVec<T> &x, NRSMat<T> &h) const
+{
+const int nq=x.size()/2;
+h.copyonwrite();
+h=(T)0;
+for(int i=0; i<nq; ++i)
+	{
+	blochqubit<T> 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<nq-1) eff+= J(i)*std::cos(x[2*i+2])*std::cos(x[2*i+3]);
+	//z_aa = z_bb = -z, x_aa = x_bb = -x, z_ab = sin a sin b, x_ab = -sin a cos b
+	h(2*i,2*i)= -eff*q.z - kfield(i)*q.x;
+	h(2*i+1,2*i+1)= h(2*i,2*i);
+	h(2*i,2*i+1)= eff*std::sin(x[2*i])*std::sin(x[2*i+1]) - kfield(i)*std::sin(x[2*i])*std::cos(x[2*i+1]);
+	if(i<nq-1)
+		{
+		blochqubit<T> 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<typename T>
+void testproblem_vqe<T>::metric(const NRVec<T> &x, NRSMat<T> &g) const
+{
+const int nq=x.size()/2;
+g.copyonwrite();
+g=(T)0;
+for(int i=0; i<nq; ++i)
+	{
+	T ca=std::cos(x[2*i]);
+	g(2*i,2*i)=0.25;
+	g(2*i+1,2*i+1)=0.25*ca*ca;
+	}
+}
+
+
 //safeguarded cubic interpolation of the minimizer of phi between a and b, given function values and derivatives at both ends
 template<typename T>
 static T cubic_interpolate(const T a, const T fa, const T da, const T b, const T fb, const T db)
@@ -428,6 +575,93 @@ return linesearch_minimize(names[hp.my_update],dir,o,x,hp,achieved);
 }
 
 
+template<typename T>
+optresult<T> GDIIS_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+const char *name="GDIIS";
+GDIIS_hyperparameters<T> hp;
+if(par) hp.setfrom(*par);
+if(hp.my_saddle_order) laerror("GDIIS cannot search for saddle points");
+if(hp.my_diis_dim<1) laerror("GDIIS needs diis_dim >= 1");
+const int verbose=hp.my_verbose;
+optevaluator<T> e(o,hp.my_maxfunc);
+optstate<T> st(name,hp,e);
+const int n=x.size();
+QN_direction<T> dir(hp,n);
+DIIS<NRVec<T>,NRVec<T> > diis(hp.my_diis_dim);
+bool fresh=true; //first step after a (re)start, the quasi-Newton step from the initial Hessian guess has no scale
+NRVec<T> g(n);
+T f=e(x,&g);
+
+for(;; ++st.iter)
+	{
+	if(st.stop(f,g.norm())) break;
+
+	NRVec<T> p = dir.direction(g); //quasi-Newton step, the error vector of DIIS
+	if(p*g >= 0)
+		{
+		if(verbose) std::cout <<name<<": not a descent direction, restart\n";
+		dir.reset();
+		diis.reset();
+		fresh=true;
+		p = g*(T)-1;
+		}
+
+	NRVec<T> xold(x), gold(g);
+	T fold=f;
+	NRVec<T> s;
+	bool accepted=false;
+	if(!fresh) //DIIS extrapolation, otherwise a line search is done directly
+		{
+		s = x+p;
+		NRVec<T> err(p);
+		diis.extrapolate(s,err,verbose>1,hp.my_diis_eps);
+		s -= x;
+		if(s*g >= 0 || s.norm() > hp.my_maxsteplength)
+			{
+			if(verbose) std::cout <<name<<": unreasonable DIIS step, restart\n";
+			diis.reset();
+			s = p;
+			T snorm=s.norm();
+			if(snorm > hp.my_maxsteplength) s *= hp.my_maxsteplength/snorm;
+			}
+		NRVec<T> xnew = x+s;
+		NRVec<T> gnew(n);
+		T fnew=e(xnew,&gnew);
+		if(fnew<f) {x=xnew; f=fnew; g=gnew; accepted=true;}
+		else
+			{
+			if(verbose) std::cout <<name<<": no decrease, restart and line search\n";
+			diis.reset();
+			}
+		}
+	if(!accepted) //line search along the quasi-Newton step
+		{
+		if(e.exhausted()) continue; //reported by stop()
+		T alpha=linesearch_wolfe(e,x,f,g,p,dir.steepest() ? 1/p.norm() : (T)1,hp);
+		if(alpha==0)
+			{
+			if(e.exhausted()) continue;
+			if(dir.steepest()) {if(verbose) std::cout <<name<<": line search failed along steepest descent, giving up\n"; break;}
+			if(verbose) std::cout <<name<<": line search failed, restart\n";
+			dir.reset();
+			fresh=true;
+			continue;
+			}
+		s = x-xold;
+		}
+	fresh=false;
+
+	NRVec<T> y = g-gold;
+	if(verbose>1) std::cout <<name<<" step |s|= "<<s.norm()<<" s.y= "<<s*y<<std::endl;
+	if(!dir.update(s,y,p,g) && verbose) std::cout <<name<<": update skipped\n";
+	st.step(f-fold,s.norm());
+	}
+
+return st.finish(achieved);
+}
+
+
 //L-BFGS: last m pairs s = x_{k+1}-x_k, y = g_{k+1}-g_k in a ring buffer
 template<typename T>
 class LBFGS_direction : public linesearch_direction<T> {
@@ -742,11 +976,42 @@ return st.finish(achieved);
 
 
 
+//natural gradient (G + reg I)^-1 g from the gradient of an inner estimator and the metric G of the problem
+template<typename T>
+class natural_gradient_estimator : public gradient_estimator<T> {
+public:
+gradient_estimator<T> &inner;
+const OptMetric<T> *m;
+T reg;
+natural_gradient_estimator(gradient_estimator<T> &in, const OptMetric<T> *mm, const T r) : inner(in), m(mm), reg(r) {};
+T operator()(optevaluator<T> &e, const NRVec<T> &x, NRVec<T> &g, const int k)
+	{
+	T f=inner(e,x,g,k);
+	NRSMat<T> G(x.size());
+	m->metric(x,G);
+	if(reg!=0) G+=reg;
+	linear_solve(G,g);
+	return f;
+	}
+int cost() const {return inner.cost();}
+bool exact() const {return inner.exact();}
+T calibrate(optevaluator<T> &e, const NRVec<T> &x, const firstorder_hyperparameters<T> &hp) {return inner.calibrate(e,x,hp);}
+};
+
+
 template<typename T>
 optresult<T> firstorder_minimize(const char *name, firstorder_step<T> &stepper, const OptProblem<T> *o, NRVec<T> &x, const firstorder_hyperparameters<T> &hp, hyperparameters<T> *achieved, gradient_estimator<T> *est)
 {
 gradient_estimator<T> exact;
-gradient_estimator<T> &grad = est ? *est : exact;
+gradient_estimator<T> *plain = est ? est : &exact;
+const OptMetric<T> *metric=NULL;
+if(hp.my_natural)
+	{
+	metric=optcapability<OptMetric<T> >(o);
+	if(!metric) laerror("natural gradient needs a problem providing the metric (OptMetric)");
+	}
+natural_gradient_estimator<T> natural(*plain,metric,hp.my_metric_reg);
+gradient_estimator<T> &grad = hp.my_natural ? natural : *plain;
 optevaluator<T> e(o,hp.my_maxfunc);
 optstate<T> st(name,hp,e);
 const int n=x.size();
@@ -802,6 +1067,142 @@ NRVec<T> step(const NRVec<T> &g, const int k)
 };
 
 
+template<typename T>
+optresult<T> QNG_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+QNG_hyperparameters<T> hp;
+if(par) hp.setfrom(*par);
+hp.my_natural=true; //also if the caller passed e.g. SGD_hyperparameters
+if(hp.my_saddle_order) laerror("QNG cannot search for saddle points");
+SGD_stepper<T> stepper(hp);
+return firstorder_minimize("QNG",stepper,o,x,hp,achieved);
+}
+
+
+//trigonometric polynomial a_0 + sum_{k=1}^R (a_k cos(k w t) + b_k sin(k w t)) of the single variable t = x[0]
+template<typename T>
+class trigpolynomial : public OptProblem<T> {
+public:
+NRVec<T> a, b;
+T w;
+trigpolynomial(const int R, const T ww) : a(R+1), b(R+1), w(ww) {};
+T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const
+	{
+	if(grad) laerror("trigpolynomial provides only the function");
+	T f=a[0];
+	for(int k=1; k<a.size(); ++k) f+= a[k]*std::cos(k*w*x[0]) + b[k]*std::sin(k*w*x[0]);
+	return f;
+	}
+};
+
+//exact minimization of f along x_i, a trigonometric polynomial of degree R with the given period: reconstruction from
+//2R+1 equidistant samples t_j = j*period/(2R+1) (the one at t=0 is the current f, 2R new evaluations) by discrete Fourier transform
+//and minimization of the reconstruction; returns false if the evaluation budget is insufficient (x,f then unchanged)
+template<typename T>
+static bool trig_minimize(optevaluator<T> &e, NRVec<T> &x, T &f, const int i, const int R, const T period, const bool reevaluate, const int verbose)
+{
+const int N=2*R+1;
+if(e.remaining()<N-1+(reevaluate?1:0)) return false;
+const T w=2*M_PI/period;
+const T xi=((const NRVec<T> &)x)[i]; //const access, x may be shared with the evaluator
+NRVec<T> fs(N); //samples
+fs[0]=f;
+NRVec<T> xs(x);
+for(int j=1; j<N; ++j)
+	{
+	xs.copyonwrite(); //the evaluator may keep a reference to the previous sample as its best point
+	xs[i]=xi+j*period/N;
+	fs[j]=e(xs);
+	}
+
+trigpolynomial<T> p(R,w);
+p.a[0]=fs.sum()/N;
+for(int k=1; k<=R; ++k)
+	{
+	p.a[k]=p.b[k]=0;
+	for(int j=0; j<N; ++j) {T phi=2*M_PI*k*j/N; p.a[k]+= fs[j]*std::cos(phi); p.b[k]+= fs[j]*std::sin(phi);}
+	p.a[k]*= (T)2/N; p.b[k]*= (T)2/N;
+	}
+
+T t, fmin;
+if(R==1) //closed form: a_1 cos(wt) + b_1 sin(wt) = A cos(wt - phi), minimum at wt = phi + pi
+	{
+	t=(std::atan2(p.b[1],p.a[1])+M_PI)/w;
+	fmin=p.a[0]-std::sqrt(p.a[1]*p.a[1]+p.b[1]*p.b[1]);
+	}
+else //global minimum of the reconstruction: best of a dense grid, refined by Brent's method (no evaluations of the problem)
+	{
+	const int M=8*N;
+	NRVec<T> tt(1);
+	int jbest=0;
+	T fbest=0;
+	for(int j=0; j<M; ++j) {tt[0]=j*period/M; T ft=p.funcgrad(tt,NULL); if(j==0 || ft<fbest) {fbest=ft; jbest=j;}}
+	optevaluator<T> pe(&p);
+	tt[0]=jbest*period/M;
+	NRVec<T> u((T)1,1);
+	linemin_brent(pe,tt,fbest,u,period/M,(T)1e-10,(T)1e-12*period,100);
+	t=tt[0];
+	fmin=fbest;
+	}
+t-= period*std::floor(t/period+0.5); //shortest step, t in (-period/2,period/2]
+x.copyonwrite();
+x[i]=xi+t;
+f= reevaluate ? e(x) : fmin;
+if(verbose>1) std::cout <<"trigonometric minimization of parameter "<<i<<" degree "<<R<<" step "<<t<<" f= "<<f<<std::endl;
+return true;
+}
+
+
+template<typename T>
+optresult<T> Rotosolve_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+const char *name="Rotosolve";
+Rotosolve_hyperparameters<T> hp;
+if(par) hp.setfrom(*par);
+if(hp.my_saddle_order) laerror("Rotosolve cannot search for saddle points");
+const int verbose=hp.my_verbose;
+const T abstol=0.01*hp.my_x_thr; //as in Powell's method
+const int n=x.size();
+const OptTrigonometric<T> *trig=optcapability<OptTrigonometric<T> >(o);
+if(!trig && verbose) std::cout <<name<<": the problem is not trigonometric, line minimizations only\n";
+optevaluator<T> e(o,hp.my_maxfunc);
+optstate<T> st(name,hp,e);
+T f=e(x);
+NRVec<T> tstep(hp.my_initstep,n); //trial steps of the line minimizations
+
+for(;; ++st.iter)
+	{
+	if(st.stop(f)) break;
+	NRVec<T> x0(x);
+	T f0=f;
+	bool budget=true;
+	for(int i=0; i<n && budget; ++i)
+		{
+		const int R= trig ? trig->degree(i) : 0;
+		if(R>0) budget=trig_minimize(e,x,f,i,R,trig->period(i),hp.my_reevaluate,verbose);
+		else
+			{
+			NRVec<T> u((T)0,n);
+			u[i]=1;
+			T t=linemin_brent(e,x,f,u,tstep[i],hp.my_linmin_tol,abstol,hp.my_maxlinemin,verbose);
+			if(t!=0) tstep[i]=std::abs(t);
+			budget=!e.exhausted();
+			}
+		}
+	NRVec<T> dx = x-x0;
+	st.step(f0-f,dx.norm());
+	if(!budget) //incomplete sweep, no convergence test
+		{
+		if(verbose) std::cout <<name<<": maximum number of function evaluations reached\n";
+		st.f=f;
+		break;
+		}
+	}
+
+return st.finish(achieved);
+}
+
+
 template<typename T>
 optresult<T> SGD_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
 {
@@ -949,12 +1350,18 @@ return firstorder_minimize("SPSA",stepper,o,x,hp,achieved,&est);
 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 testproblem_vqe<double>;
+template class ParameterShift<double>;
+template class FiniteDifference<double>;
 template class BFGS_Optimizer<double>;
 template class LBFGS_Optimizer<double>;
+template class GDIIS_Optimizer<double>;
 template class CG_Optimizer<double>;
 template class NelderMead_Optimizer<double>;
 template class Powell_Optimizer<double>;
+template class Rotosolve_Optimizer<double>;
 template class SGD_Optimizer<double>;
+template class QNG_Optimizer<double>;
 template class Adam_Optimizer<double>;
 template class SPSA_Optimizer<double>;
 template double linemin_brent(optevaluator<double> &e, NRVec<double> &x, double &f, const NRVec<double> &u, const double t1, const double tol, const double abstol, const int maxiter, const int verbose);
diff --git a/optimize.h b/optimize.h
index cfd6d3f..15892bf 100644
--- a/optimize.h
+++ b/optimize.h
@@ -24,6 +24,7 @@
 
 #include <stdint.h>
 #include <cstdarg>
+#include <cmath>
 #include "vec.h"
 #include "mat.h"
 #include "smat.h"
@@ -41,6 +42,79 @@ virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const =0; //evaluate funct
 virtual ~OptProblem() {};
 };
 
+//optional capabilities of a problem, to be multiply inherited together with OptProblem; methods find them by optcapability()
+//matrices are allocated by the caller, implementations writing their elements have to call copyonwrite() first
+
+//Hessian of f
+template<typename T>
+class OptHessian {
+public:
+virtual void hessian(const NRVec<T> &x, NRSMat<T> &h) const =0; //Hessian at x
+virtual NRVec<T> hessian_times(const NRVec<T> &x, const NRVec<T> &v) const {NRSMat<T> h(x.size()); hessian(x,h); return h*v;} //Hessian-vector product, can be overridden for large problems
+virtual ~OptHessian() {};
+};
+
+//metric tensor of the parameter space: Fubini-Study metric for the quantum natural gradient (Stokes et al., Quantum 4, 269 (2020)),
+//Fisher information for the natural gradient in machine learning (Amari, Neural Comput. 10, 251 (1998))
+template<typename T>
+class OptMetric {
+public:
+virtual void metric(const NRVec<T> &x, NRSMat<T> &g) const =0; //metric at x
+virtual ~OptMetric() {};
+};
+
+//f as a function of x_i alone is a trigonometric polynomial a_0 + sum_{k=1}^{R_i} (a_k cos(k w_i x_i) + b_k sin(k w_i x_i)), w_i = 2 pi/period_i;
+//for parameterized quantum circuits with gates exp(-i x_i P/2), P a Pauli string, R_i = 1 and period_i = 2 pi
+template<typename T>
+class OptTrigonometric {
+public:
+virtual int degree(const int i) const =0; //R_i, 0 = f is not trigonometric in x_i
+virtual T period(const int i) const {return 2*M_PI;}
+virtual ~OptTrigonometric() {};
+};
+
+
+//problem wrapping another one (e.g. providing gradients); the capabilities of the wrapped problem remain accessible by optcapability()
+template<typename T>
+class OptProblemAdapter : public OptProblem<T> {
+public:
+const OptProblem<T> *p; //wrapped problem
+OptProblemAdapter(const OptProblem<T> *pp) : p(pp) {};
+};
+
+//capability C (OptHessian<T>, OptMetric<T>, ...) of the problem, searched also through adapters; NULL if not available
+template<class C, typename T>
+const C *optcapability(const OptProblem<T> *o)
+{
+const C *c=dynamic_cast<const C *>(o);
+if(c) return c;
+const OptProblemAdapter<T> *a=dynamic_cast<const OptProblemAdapter<T> *>(o);
+return a ? optcapability<C>(a->p) : NULL;
+}
+
+
+//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
+template<typename T>
+class ParameterShift : public OptProblemAdapter<T> {
+public:
+const OptTrigonometric<T> *trig;
+T h;
+ParameterShift(const OptProblem<T> *pp, const T hh=1e-6);
+virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const;
+};
+
+//gradient by central differences with step h, 2 evaluations per parameter
+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) {};
+virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const;
+};
+
+
 //test case
 template<typename T>
 class testproblem : public OptProblem<T> {
@@ -49,6 +123,22 @@ virtual T funcgrad(const NRVec<T> &x, NRVec<T> *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<typename T>
+class testproblem_vqe : public OptProblem<T>, public OptHessian<T>, public OptMetric<T>, public OptTrigonometric<T> {
+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<T> &x, NRVec<T> *grad) const;
+virtual void hessian(const NRVec<T> &x, NRSMat<T> &h) const;
+virtual void metric(const NRVec<T> &x, NRSMat<T> &g) const;
+virtual int degree(const int i) const {return 1;}
+virtual ~testproblem_vqe() {};
+};
+
 
 //base classes for optimization methods
 
@@ -206,7 +296,9 @@ T my_lr_decay; //learning rate in iteration k is lr/(1+lr_decay*k)^lr_power
 T my_lr_power;
 T my_weight_decay; //decoupled weight decay x -= lr*weight_decay*x (Loshchilov & Hutter, ICLR 2019)
 bool my_return_best; //return the point with the lowest function value found instead of the last one; only for (nearly) deterministic f, the lowest of noisy values is biased by the noise
-firstorder_hyperparameters() : my_lr(1e-3), my_lr_decay(0), my_lr_power(1), my_weight_decay(0), my_return_best(false) {this->my_maxstep=1000;};
+bool my_natural; //natural gradient (G + metric_reg I)^-1 g with the metric G of the problem (OptMetric) (Amari, Neural Comput. 10, 251 (1998); Stokes et al., Quantum 4, 269 (2020)); |g| refers then to it
+T my_metric_reg; //regularization added to the diagonal of the metric
+firstorder_hyperparameters() : my_lr(1e-3), my_lr_decay(0), my_lr_power(1), my_weight_decay(0), my_return_best(false), my_natural(false), my_metric_reg(1e-3) {this->my_maxstep=1000;};
 virtual void setfrom(const hyperparameters<T> &p) {const firstorder_hyperparameters *q=dynamic_cast<const firstorder_hyperparameters *>(&p); if(q) *this= *q; else gradient_hyperparameters<T>::setfrom(p);}
 };
 
@@ -284,6 +376,27 @@ virtual ~LBFGS_Optimizer() {};
 };
 
 
+//geometry DIIS (Csaszar & Pulay, J. Mol. Struct. 114, 31 (1984)) with a quasi-Newton Hessian (Farkas & Schlegel, Phys. Chem. Chem. Phys. 4, 11 (2002))
+//error vectors are the quasi-Newton steps e_i = -H g_i, the new point is sum_i c_i (x_i+e_i) minimizing |sum_i c_i e_i| with sum_i c_i = 1
+//one function evaluation per iteration; if the DIIS step goes uphill, is longer than maxsteplength or does not decrease f,
+//the DIIS history is restarted and a line search along the quasi-Newton step is done
+template<typename T>
+class GDIIS_hyperparameters : public BFGS_hyperparameters<T> {
+public:
+int my_diis_dim; //no. of stored points
+T my_diis_eps; //regularization added to the diagonal of the DIIS matrix
+GDIIS_hyperparameters() : my_diis_dim(5), my_diis_eps(0) {};
+virtual void setfrom(const hyperparameters<T> &p) {const GDIIS_hyperparameters *q=dynamic_cast<const GDIIS_hyperparameters *>(&p); if(q) *this= *q; else BFGS_hyperparameters<T>::setfrom(p);}
+};
+
+template<typename T>
+class GDIIS_Optimizer : public Optimizer<T> {
+public:
+virtual optresult<T> minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par=NULL, hyperparameters<T> *achieved=NULL);
+virtual ~GDIIS_Optimizer() {};
+};
+
+
 //nonlinear conjugate gradient methods (Nocedal & Wright chap. 5.2; Hager & Zhang, SIAM J. Optim. 16, 170 (2005))
 //differing by the formula for beta in p_{k+1} = -g_{k+1} + beta p_k; PR and HS are used in the nonnegative variants PR+ and HS+
 enum CG_beta {CG_FR, CG_PR, CG_HS, CG_DY, CG_HZ};
@@ -349,6 +462,26 @@ virtual ~Powell_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;
+//parameters of degree 0, or all if the problem is not trigonometric, are minimized by line minimization as in Powell's method (whose parameters are used)
+template<typename T>
+class Rotosolve_hyperparameters : public Powell_hyperparameters<T> {
+public:
+bool my_reevaluate; //evaluate f at the new point instead of taking the value of the reconstruction (for noisy f or doubtful degrees)
+Rotosolve_hyperparameters() : my_reevaluate(false) {};
+virtual void setfrom(const hyperparameters<T> &p) {const Rotosolve_hyperparameters *q=dynamic_cast<const Rotosolve_hyperparameters *>(&p); if(q) *this= *q; else Powell_hyperparameters<T>::setfrom(p);}
+};
+
+template<typename T>
+class Rotosolve_Optimizer : public Optimizer<T> {
+public:
+virtual optresult<T> minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par=NULL, hyperparameters<T> *achieved=NULL);
+virtual ~Rotosolve_Optimizer() {};
+};
+
+
 //stochastic gradient descent with momentum v = momentum*v + g, step d = v or d = g + momentum*v (Nesterov) (Sutskever et al., ICML 2013)
 template<typename T>
 class SGD_hyperparameters : public firstorder_hyperparameters<T> {
@@ -367,6 +500,24 @@ virtual ~SGD_Optimizer() {};
 };
 
 
+//quantum natural gradient x -= lr (G + metric_reg I)^-1 g (Stokes et al., Quantum 4, 269 (2020)): SGD without momentum with the natural gradient;
+//the problem has to provide the metric (OptMetric), the gradient can come e.g. from ParameterShift; momentum, Adam or SPSA can be combined with the natural gradient by setting my_natural
+//natural steps are longer than gradient steps with the same lr (metric eigenvalues ~1/4 or smaller): lr up to ~0.1, or a larger metric_reg
+template<typename T>
+class QNG_hyperparameters : public SGD_hyperparameters<T> {
+public:
+QNG_hyperparameters() {this->my_momentum=0; this->my_natural=true; this->my_lr=0.05;};
+virtual void setfrom(const hyperparameters<T> &p) {const QNG_hyperparameters *q=dynamic_cast<const QNG_hyperparameters *>(&p); if(q) *this= *q; else SGD_hyperparameters<T>::setfrom(p);}
+};
+
+template<typename T>
+class QNG_Optimizer : public Optimizer<T> {
+public:
+virtual optresult<T> minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par=NULL, hyperparameters<T> *achieved=NULL);
+virtual ~QNG_Optimizer() {};
+};
+
+
 //Adam (Kingma & Ba, ICLR 2015), with weight_decay>0 AdamW, optionally AMSGrad (Reddi et al., ICLR 2018); beta1=0 gives RMSprop
 template<typename T>
 class Adam_hyperparameters : public firstorder_hyperparameters<T> {
diff --git a/t.cc b/t.cc
index e5585f0..fe62ec2 100644
--- a/t.cc
+++ b/t.cc
@@ -4830,7 +4830,7 @@ cout <<"adjoint error = "<<(a*adj).norm(abs(det))<<endl;
 
 }
 
-if(1)
+if(0) //testbed of the general methods
 {
 int seed;
 int f=open("/dev/random",O_RDONLY);
@@ -4850,12 +4850,13 @@ testproblem<double> problem;
 //Powell_hyperparameters<double> param; //tested
 //SGD_hyperparameters<double> param; //tested
 //Adam_hyperparameters<double> param; //tested
-SPSA_hyperparameters<double> param;
+//SPSA_hyperparameters<double> param; //tested
+GDIIS_hyperparameters<double> param;
 param.my_maxstep=100;
 param.my_f_thr=1e-6;
 param.my_x_thr=1e-3;
 param.my_g_thr=1e-3; //gradient methods
-param.my_lr=1e-2; //first-order methods
+//param.my_lr=1e-2; //first-order methods
 param.my_maxstep=2000; //first-order methods
 //param.my_update=QN_BFGS; //tested
 //param.my_update=QN_DFP; //tested
@@ -4868,7 +4869,8 @@ cout << "initial " <<x<<endl;
 //Powell_Optimizer<double> o; //tested
 //SGD_Optimizer<double> o; //tested
 //Adam_Optimizer<double> o; //tested
-SPSA_Optimizer<double> o;
+//SPSA_Optimizer<double> o; //tested
+GDIIS_Optimizer<double> o;
 gradient_hyperparameters<double> achieved; //gradient methods
 //hyperparameters<double> achieved; //derivative-free methods
 optresult<double> res=o.minimize(&problem,x,&param,&achieved);
@@ -4879,4 +4881,38 @@ double r=problem.funcgrad(x,NULL);
 cout <<r<<endl;
 }
 
+if(1) //testbed of the VQE methods, the problem provides capabilities (OptMetric, OptTrigonometric, OptHessian)
+{
+int seed;
+int f=open("/dev/random",O_RDONLY);
+if(sizeof(int)!=read(f,&seed,sizeof(int))) laerror("cannot read /dev/random");
+close(f);
+srand(seed);
+
+int n;
+cin >>n;
+NRVec<double> x(n);
+x.randomize(1.);
+testproblem_vqe<double> problem; //needs an even dimension
+ParameterShift<double> psproblem(&problem); //gradient by the parameter-shift rule from function values only
+//QNG_hyperparameters<double> param; //tested
+Rotosolve_hyperparameters<double> 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 " <<x<<endl;
+//QNG_Optimizer<double> o; //tested
+Rotosolve_Optimizer<double> o;
+//gradient_hyperparameters<double> achieved; //gradient methods
+hyperparameters<double> achieved; //derivative-free methods
+optresult<double> res=o.minimize(&problem,x,&param,&achieved);
+//optresult<double> res=o.minimize(&psproblem,x,&param,&achieved); //gradient methods with parameter-shift gradients
+cout << "optimized " <<x<<endl;
+//cout << "converged "<<res.converged<<" f= "<<res.f<<" iterations "<<achieved.my_maxstep<<" function evaluations "<<achieved.my_maxfunc<<" last df "<<achieved.my_f_thr<<" last dx "<<achieved.my_x_thr<<" |g| "<<achieved.my_g_thr<<endl; //gradient methods
+cout << "converged "<<res.converged<<" f= "<<res.f<<" iterations "<<achieved.my_maxstep<<" function evaluations "<<achieved.my_maxfunc<<" last df "<<achieved.my_f_thr<<" last dx "<<achieved.my_x_thr<<endl; //derivative-free methods
+double r=problem.funcgrad(x,NULL);
+cout <<r<<endl;
+}
+
 }//main