diff --git a/optimize.cc b/optimize.cc
index 8c300e9..f0c843e 100644
--- a/optimize.cc
+++ b/optimize.cc
@@ -83,6 +83,47 @@ template<typename T>
 void FrozenOptProblem<T>::constraints(const NRVec<T> &x, NRVec<T> &c) const {optcapability<OptConstraints<T> >(this->p)->constraints(expand(x),c);}
 
 
+//Jacobian of the residuals by forward differences, columns (r(x + h_j e_j) - r(x))/h_j, h_j = h max(1,|x_j|); resid(x,r) evaluates the residuals
+template<typename T, typename F>
+static void forward_difference_jacobian(F resid, const NRVec<T> &x, const NRVec<T> &r, NRMat<T> &J, const T h)
+{
+NRVec<T> xs(x);
+for(int j=0; j<x.size(); ++j)
+	{
+	const T hj=h*std::max((T)1,std::abs(x[j]));
+	xs.copyonwrite(); //the evaluator may keep a reference to the previous point as its best one
+	xs[j]=x[j]+hj;
+	NRVec<T> rs(r.size());
+	resid(xs,rs);
+	rs-=r;
+	rs*= 1/hj;
+	J.columnset(rs,j);
+	xs.copyonwrite();
+	xs[j]=x[j];
+	}
+}
+
+template<typename T>
+void OptLeastSquares<T>::jacobian_fd(const NRVec<T> &x, const NRVec<T> &r, NRMat<T> &J) const
+{
+forward_difference_jacobian([this](const NRVec<T> &xx, NRVec<T> &rr){this->residuals(xx,rr);},x,r,J,fd_step);
+}
+
+template<typename T>
+T OptLeastSquares<T>::funcgrad(const NRVec<T> &x, NRVec<T> *grad) const
+{
+NRVec<T> r(this->nresiduals());
+this->residuals(x,r);
+if(grad)
+	{
+	NRMat<T> J(r.size(),x.size());
+	if(!this->jacobian(x,J)) jacobian_fd(x,r,J);
+	*grad = r*J; //J^T r
+	}
+return (r*r)/2;
+}
+
+
 //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)
@@ -734,6 +775,402 @@ return linesearch_minimize(names[hp.my_beta],dir,o,x,hp,achieved);
 
 
 
+template<typename T>
+bool hessian_update(NRSMat<T> &B, const NRVec<T> &s, const NRVec<T> &y, const TR_hessian type, const T thr)
+{
+const NRVec<T> Bs=B*s;
+const T ss=s*s;
+NRVec<T> xi(y);
+xi-=Bs; //y - Bs, vanishes if B already satisfies the secant condition
+const T xis=xi*s;
+switch(type)
+	{
+	case TR_BFGS: //B+ = B - Bs (Bs)^T/(s.Bs) + y y^T/(y.s)
+		{
+		const T sy=s*y, sBs=s*Bs;
+		if(sy<=thr*std::sqrt(ss*(y*y)) || sBs<=0) return false;
+		B.rank2update(1/(2*sy),y,y);
+		B.rank2update(-1/(2*sBs),Bs,Bs);
+		}
+		break;
+	case TR_SR1: //B+ = B + xi xi^T/(xi.s) (Nocedal & Wright eq. 6.24)
+		if(std::abs(xis)<=thr*xi.norm()*std::sqrt(ss)) return false;
+		B.rank2update(1/(2*xis),xi,xi);
+		break;
+	case TR_PSB: //B+ = B + (xi s^T + s xi^T)/(s.s) - (xi.s) s s^T/(s.s)^2 (Powell-symmetric-Broyden)
+		B.rank2update(1/ss,xi,s);
+		B.rank2update(-xis/(2*ss*ss),s,s);
+		break;
+	case TR_BOFILL: //phi*SR1 + (1-phi)*PSB, phi = (xi.s)^2/(|xi|^2 |s|^2), the SR1 term phi xi xi^T/(xi.s) = (xi.s) xi xi^T/(|xi|^2 |s|^2) has no small denominator
+		{
+		const T xixi=xi*xi;
+		if(xixi==0) return false;
+		const T phi=xis*xis/(xixi*ss);
+		B.rank2update(xis/(2*xixi*ss),xi,xi);
+		B.rank2update((1-phi)/ss,xi,s);
+		B.rank2update(-(1-phi)*xis/(2*ss*ss),s,s);
+		}
+		break;
+	default: laerror("unknown Hessian update");
+	}
+return true;
+}
+
+
+template<typename T>
+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)
+{
+const int verbose=hp.my_verbose;
+const bool saddle=stepper.saddle();
+const OptHessian<T> *hes=optcapability<OptHessian<T> >(o);
+const bool exact= hp.my_hessian==TR_EXACT || (hp.my_hessian==TR_AUTO && hes);
+if(exact && !hes) laerror("exact Hessian requested, but the problem does not provide OptHessian");
+TR_hessian update=hp.my_hessian; //update of B between exact Hessians or without them: Bofill for saddle points; for minima SR1 after
+if(update==TR_AUTO || update==TR_EXACT) update= saddle ? TR_BOFILL : (exact ? TR_SR1 : TR_BFGS); //an exact Hessian, which may be indefinite, BFGS from the identity
+if(saddle && update==TR_BFGS) laerror("the BFGS update keeps the Hessian positive definite, unsuitable for saddle points");
+optevaluator<T> e(o,hp.my_maxfunc);
+optstate<T> st(name,hp,e);
+const int n=x.size();
+NRVec<T> g(n);
+T f=e(x,&g);
+NRSMat<T> B(n);
+bool Binit=false; //B is the unscaled initial guess
+if(exact) hes->hessian(x,B);
+else {B=1.; Binit=true;}
+int sinceexact=0;
+T delta=hp.my_delta0;
+
+for(;; ++st.iter)
+	{
+	if(st.stop(f,g.norm())) break;
+	NRVec<T> p=stepper.step(g,B,delta);
+	const T pnorm=p.norm();
+	NRVec<T> Bp=B*p;
+	const T pred= g*p + 0.5*(p*Bp); //predicted change of f
+	NRVec<T> xn=x+p;
+	NRVec<T> gn(n);
+	const T fn=e(xn,&gn);
+	const T ratio= pred!=0 ? (fn-f)/pred : 0;
+	if(verbose>1) std::cout <<name<<" step |p|= "<<pnorm<<" delta= "<<delta<<" f= "<<fn<<" ratio= "<<ratio<<std::endl;
+	bool accept;
+	if(saddle) //the step is always taken, the radius follows the quality of the model
+		{
+		accept=true;
+		if(ratio<0.25 || ratio>1.75) delta=0.25*pnorm;
+		else if(ratio>0.75 && ratio<1.25 && pnorm>=0.99*delta) delta=std::min(2*delta,hp.my_deltamax);
+		}
+	else if(std::abs(pred)<=100*std::numeric_limits<T>::epsilon()*std::abs(f)) //change of f below its rounding errors, judged by the gradient
+		{
+		accept= gn.norm()<g.norm();
+		if(!accept) delta=0.25*pnorm;
+		}
+	else //Nocedal & Wright alg. 4.1
+		{
+		accept= pred<0 && ratio>hp.my_eta;
+		if(ratio<0.25) delta=0.25*pnorm;
+		else if(ratio>0.75 && pnorm>=0.99*delta) delta=std::min(2*delta,hp.my_deltamax);
+		}
+	if(!accept)
+		{
+		if(verbose) std::cout <<name<<": step rejected\n";
+		continue;
+		}
+	NRVec<T> y=gn-g;
+	const T fold=f;
+	x=xn; f=fn; g=gn;
+	if(exact && ++sinceexact>=hp.my_hessian_recompute) {hes->hessian(x,B); sinceexact=0;}
+	else
+		{
+		if(Binit) //scaling of the initial guess (cf. Nocedal & Wright eq. 6.20)
+			{
+			const T sy=p*y;
+			if(sy>0) B=(T)((y*y)/sy);
+			Binit=false;
+			}
+		if(!hessian_update(B,p,y,update,hp.my_update_thr) && verbose>1) std::cout <<name<<": Hessian update skipped\n";
+		}
+	st.step(f-fold,pnorm);
+	}
+
+//number of negative eigenvalues of the Hessian at the result
+NRSMat<T> H(n);
+if(hes) hes->hessian(x,H);
+else H=B;
+NRVec<T> w(n);
+diagonalize(H,w,(NRMat<T> *)NULL);
+int index=0;
+for(int i=0; i<n; ++i) if(w[i]<0) ++index;
+if(verbose) std::cout <<name<<": "<<index<<" negative eigenvalues of the"<<(hes ? "" : " approximate")<<" Hessian\n";
+optresult<T> r=st.finish(achieved);
+if(achieved) achieved->my_saddle_order=index;
+return r;
+}
+
+
+//-a/(lambda+mu), for the components of a trust-region step in the eigenbasis of B
+template<typename T>
+class shifteddivision {
+public:
+T mu;
+shifteddivision(const T m) : mu(m) {};
+T operator()(const T &a, const T &lambda) const {return -a/(lambda+mu);}
+};
+
+//exact solution of min g.p + 1/2 p.B.p, |p| <= delta in the eigenbasis of B = V Lambda V^T: p(mu) = -sum_i (v_i.g)/(lambda_i+mu) v_i
+//with mu = 0 for the Newton step inside the trust region, otherwise mu > -lambda_min from |p(mu)| = delta (More & Sorensen 1983)
+template<typename T>
+static NRVec<T> trustregion_eigenstep(const NRVec<T> &g, const NRSMat<T> &B, const T delta)
+{
+const int n=g.size();
+NRSMat<T> Bc(B);
+NRVec<T> lam(n);
+NRMat<T> V(n,n);
+diagonalize(Bc,lam,&V); //ascending eigenvalues, eigenvectors in the columns of V
+const NRVec<T> gt=g*V; //components of g in the eigenbasis
+const T lmin=lam[0];
+if(lmin>0) //Newton step if inside the trust region
+	{
+	NRVec<T> pt(gt);
+	pt.call_on_me(shifteddivision<T>(0),lam);
+	if(pt.norm()<=delta) return V*pt;
+	}
+const T lo=std::max((T)0,-lmin);
+const T tol=1e-12*std::max(std::abs(lam[0]),std::abs(lam[n-1]));
+//hard case: g has no component along the eigenvectors of lambda_min and the step from the others at mu=-lambda_min is inside the trust region
+if(lmin<=0)
+	{
+	int nlow=0;
+	while(nlow<n && lam[nlow]-lmin<=tol) ++nlow;
+	NRVec<T> gtlow=gt.subvector(0,nlow-1);
+	if(gtlow.norm()<=1e-12*g.norm())
+		{
+		NRVec<T> pt(gt);
+		pt.copyonwrite();
+		for(int i=0; i<nlow; ++i) pt[i]=0;
+		for(int i=nlow; i<n; ++i) pt[i]= -gt[i]/(lam[i]+lo);
+		T pn=pt.norm();
+		if(pn<=delta)
+			{
+			NRVec<T> p=V*pt;
+			p.axpy(std::sqrt(delta*delta-pn*pn),V.column(0));
+			return p;
+			}
+		}
+	}
+//|p(mu)| decreases from infinity at mu=lo; root by Newton's method on 1/|p(mu)| - 1/delta, safeguarded by bisection
+T mulo=lo, muhi=lo+g.norm()/delta; //|p(muhi)| <= delta
+T mu=muhi;
+NRVec<T> pt(gt);
+for(int it=0; it<100; ++it)
+	{
+	pt=gt;
+	pt.call_on_me(shifteddivision<T>(mu),lam);
+	const T pn=pt.norm();
+	if(std::abs(pn-delta)<=1e-10*delta) break;
+	if(pn>delta) mulo=mu; else muhi=mu;
+	NRVec<T> qt(pt); //|q|^2 = sum_i pt_i^2/(lambda_i+mu) = -d|p|^2/dmu /2
+	qt.call_on_me(shifteddivision<T>(mu),lam);
+	const T qq= -(pt*qt);
+	T munew= qq>0 ? mu+(pn*pn/qq)*(pn-delta)/delta : mu;
+	if(!(munew>mulo && munew<muhi)) munew=(mulo+muhi)/2;
+	mu=munew;
+	}
+return V*pt;
+}
+
+//approximate solution of the trust-region subproblem by truncated conjugate gradients (Steihaug 1983; Nocedal & Wright alg. 7.2)
+template<typename T>
+static NRVec<T> trustregion_steihaug(const NRVec<T> &g, const NRSMat<T> &B, const T delta, int maxit)
+{
+const int n=g.size();
+if(maxit<=0) maxit=n;
+NRVec<T> z((T)0,n);
+const T gnorm=g.norm();
+if(gnorm==0) return z;
+const T eps=std::min((T)0.5,std::sqrt(gnorm))*gnorm;
+NRVec<T> r(g);
+NRVec<T> d=g*(T)-1;
+for(int j=0; j<maxit; ++j)
+	{
+	NRVec<T> Bd=B*d;
+	const T dBd=d*Bd;
+	const T rr=r*r;
+	NRVec<T> z1(z);
+	if(dBd>0) z1.axpy(rr/dBd,d);
+	if(dBd<=0 || z1.norm()>=delta) //negative curvature or outside: to the boundary along d, |z + tau d| = delta
+		{
+		const T a=d*d, b=2*(z*d), c=z*z-delta*delta;
+		z.axpy((-b+std::sqrt(b*b-4*a*c))/(2*a),d);
+		return z;
+		}
+	r.axpy(rr/dBd,Bd);
+	z=z1;
+	if(r.norm()<eps) return z;
+	d*= (r*r)/rr;
+	d-=r;
+	}
+return z;
+}
+
+template<typename T>
+class TR_newtonstep : public trustregion_step<T> {
+public:
+const TrustRegion_hyperparameters<T> &hp;
+TR_newtonstep(const TrustRegion_hyperparameters<T> &h) : hp(h) {};
+NRVec<T> step(const NRVec<T> &g, const NRSMat<T> &B, const T delta) {return hp.my_subproblem==TR_STEIHAUG ? trustregion_steihaug(g,B,delta,hp.my_maxcg) : trustregion_eigenstep(g,B,delta);}
+};
+
+template<typename T>
+optresult<T> TrustRegion_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+TrustRegion_hyperparameters<T> hp;
+if(par) hp.setfrom(*par);
+if(hp.my_saddle_order) laerror("TrustRegion_Optimizer minimizes, use PRFO_Optimizer for saddle points");
+TR_newtonstep<T> stepper(hp);
+return trustregion_minimize(hp.my_subproblem==TR_STEIHAUG ? "TR-Steihaug" : "TR-Newton",stepper,o,x,hp,achieved);
+}
+
+
+//shift of the RFO step for the modes from..to-1: largest (maximization) or lowest (minimization) eigenvalue of the bordered matrix [[Lambda, gt],[gt^T, 0]]
+template<typename T>
+static T rfo_shift(const NRVec<T> &lam, const NRVec<T> &gt, const int from, const int to, const bool largest)
+{
+const int m=to-from;
+NRSMat<T> A(m+1);
+A=(T)0;
+for(int i=0; i<m; ++i) {A(i,i)=lam[from+i]; A(m,i)=gt[from+i];}
+NRVec<T> w(m+1);
+diagonalize(A,w,(NRMat<T> *)NULL);
+return largest ? w[m] : w[0];
+}
+
+//P-RFO step: in the eigenbasis of B, p_i = -gt_i/(lambda_i - mu) with mu from the RFO matrices of the k maximized and of the minimized modes,
+//scaled to the trust region if longer
+template<typename T>
+class PRFO_step : public trustregion_step<T> {
+public:
+int k; //no. of maximized modes
+PRFO_step(const int kk) : k(kk) {};
+bool saddle() const {return k>0;}
+NRVec<T> step(const NRVec<T> &g, const NRSMat<T> &B, const T delta)
+	{
+	const int n=g.size();
+	NRSMat<T> Bc(B);
+	NRVec<T> lam(n);
+	NRMat<T> V(n,n);
+	diagonalize(Bc,lam,&V);
+	const NRVec<T> gt=g*V;
+	const T mumax= k>0 ? rfo_shift(lam,gt,0,k,true) : 0;
+	const T mumin= k<n ? rfo_shift(lam,gt,k,n,false) : 0;
+	NRVec<T> pt(n);
+	for(int i=0; i<n; ++i)
+		{
+		const T den=lam[i]-(i<k ? mumax : mumin);
+		pt[i]= den!=0 ? -gt[i]/den : 0;
+		}
+	NRVec<T> p=V*pt;
+	const T pn=p.norm();
+	if(pn>delta) p*= delta/pn;
+	return p;
+	}
+};
+
+template<typename T>
+optresult<T> PRFO_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+PRFO_hyperparameters<T> hp;
+if(par) hp.setfrom(*par);
+if(hp.my_saddle_order<0 || hp.my_saddle_order>x.size()) laerror("invalid saddle order in PRFO_Optimizer");
+PRFO_step<T> stepper(hp.my_saddle_order);
+return trustregion_minimize("P-RFO",stepper,o,x,hp,achieved);
+}
+
+
+//Jacobian for LM: from the problem, or by forward differences with counted evaluations; false if the evaluation budget is insufficient
+template<typename T>
+static bool lm_jacobian(optevaluator<T> &e, const OptResiduals<T> *res, const NRVec<T> &x, const NRVec<T> &r, NRMat<T> &J, const T h)
+{
+if(res->jacobian(x,J)) return true;
+if(e.remaining()<x.size()) return false;
+forward_difference_jacobian([&e,res](const NRVec<T> &xx, NRVec<T> &rr){e(res,xx,rr);},x,r,J,h);
+return true;
+}
+
+template<typename T>
+optresult<T> LM_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+const char *name="Levenberg-Marquardt";
+LM_hyperparameters<T> hp;
+if(par) hp.setfrom(*par);
+if(hp.my_saddle_order) laerror("Levenberg-Marquardt cannot search for saddle points");
+const OptResiduals<T> *res=optcapability<OptResiduals<T> >(o);
+if(!res) laerror("Levenberg-Marquardt needs a least-squares problem providing OptResiduals");
+const int verbose=hp.my_verbose;
+const int n=x.size(), m=res->nresiduals();
+optevaluator<T> e(o,hp.my_maxfunc);
+optstate<T> st(name,hp,e);
+NRVec<T> r(m);
+T f=e(res,x,r);
+NRMat<T> J(m,n);
+const bool fdjac= !res->jacobian(x,J); //the problem does not provide the Jacobian
+bool budget= !fdjac || lm_jacobian(e,res,x,r,J,hp.my_fd_step);
+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;
+
+for(; budget; ++st.iter)
+	{
+	if(st.stop(f,g.norm(),fdjac ? n+1 : 1)) break;
+	//damped step (J^T J + mu D) h = -J^T r
+	NRSMat<T> M(A);
+	M.copyonwrite();
+	for(int i=0; i<n; ++i) M(i,i)+= mu*(hp.my_scaling ? diag[i] : (T)1);
+	NRVec<T> h=g*(T)-1;
+	linear_solve(M,h);
+	NRVec<T> xn=x+h;
+	NRVec<T> rn(m);
+	const T fn=e(res,xn,rn);
+	NRVec<T> Dh(h);
+	if(hp.my_scaling) Dh*=diag;
+	Dh*=mu;
+	Dh-=g;
+	const T pred=(h*Dh)/2; //predicted reduction L(0)-L(h) of the linear model
+	const T rho=(f-fn)/pred;
+	if(verbose>1) std::cout <<name<<" step |h|= "<<h.norm()<<" mu= "<<mu<<" f= "<<fn<<" gain ratio= "<<rho<<std::endl;
+	if(rho>0)
+		{
+		const T fold=f;
+		x=xn; r=rn; f=fn;
+		NRMat<T> Jn(m,n);
+		budget=lm_jacobian(e,res,x,r,Jn,hp.my_fd_step);
+		J=Jn;
+		A=J.transposedtimes();
+		g=r*J;
+		for(int i=0; i<n; ++i) diag[i]=((const NRSMat<T> &)A)(i,i);
+		mu*= std::max((T)1/3,1-std::pow(2*rho-1,3));
+		nu=2;
+		st.step(fold-f,h.norm());
+		}
+	else
+		{
+		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)
+			{
+			if(verbose) std::cout <<name<<": step at machine precision\n";
+			st.converged=true;
+			break;
+			}
+		}
+	}
+if(!budget && verbose) std::cout <<name<<": maximum number of function evaluations reached\n";
+st.f=f;
+return st.finish(achieved);
+}
+
+
 //simplex of n+1 vertices with their function values
 template<typename T>
 class NM_simplex {
@@ -1616,6 +2053,11 @@ template class FrozenOptProblem<double>;
 template class BFGS_Optimizer<double>;
 template class LBFGS_Optimizer<double>;
 template class GDIIS_Optimizer<double>;
+template class TrustRegion_Optimizer<double>;
+template class PRFO_Optimizer<double>;
+template class LM_Optimizer<double>;
+template class OptLeastSquares<double>;
+template bool hessian_update(NRSMat<double> &B, const NRVec<double> &s, const NRVec<double> &y, const TR_hessian type, const double thr);
 template class CG_Optimizer<double>;
 template class NelderMead_Optimizer<double>;
 template class Powell_Optimizer<double>;
diff --git a/optimize.h b/optimize.h
index 4f8e93d..463b8ea 100644
--- a/optimize.h
+++ b/optimize.h
@@ -85,6 +85,29 @@ virtual ~OptConstraints() {};
 };
 
 
+//least-squares problem f(x) = 1/2 |r(x)|^2 given by the residuals r and optionally their Jacobian
+template<typename T>
+class OptResiduals {
+public:
+virtual int nresiduals() const =0;
+virtual void residuals(const NRVec<T> &x, NRVec<T> &r) const =0; //r is allocated by the caller with nresiduals() elements
+virtual bool jacobian(const NRVec<T> &x, NRMat<T> &J) const {return false;} //J (nresiduals x n) allocated by the caller; false = not provided, finite differences are used
+virtual ~OptResiduals() {};
+};
+
+//base class for least-squares problems: a user derives from it and writes residuals() (and optionally jacobian());
+//funcgrad() provides f = 1/2 |r|^2 and its gradient J^T r, so that all methods can be used, not only LM_Optimizer
+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) {};
+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() {};
+};
+
+
 //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> {
@@ -229,6 +252,13 @@ T operator()(const NRVec<T> &x, NRVec<T> *grad=NULL)
 	if(nfunc++ ==0 || f<fbest) {fbest=f; xbest=x;}
 	return f;
 	}
+T operator()(const OptResiduals<T> *res, const NRVec<T> &x, NRVec<T> &r) //residuals at x, returns f = 1/2 |r|^2
+	{
+	res->residuals(x,r);
+	T f=(r*r)/2;
+	if(nfunc++ ==0 || f<fbest) {fbest=f; xbest=x;}
+	return f;
+	}
 bool exhausted() const {return maxfunc>0 && nfunc>=maxfunc;}
 int remaining() const {return maxfunc>0 ? maxfunc-nfunc : std::numeric_limits<int>::max();}
 };
@@ -433,6 +463,107 @@ virtual ~CG_Optimizer() {};
 };
 
 
+//trust-region methods: quadratic model f + g.p + 1/2 p.B.p minimized (or for saddle points stationary) within |p| <= delta
+//(Nocedal & Wright chap. 4; Conn, Gould & Toint, Trust-Region Methods, SIAM 2000)
+
+//Hessian B of the model: exact (OptHessian of the problem) or updated by BFGS, SR1, Powell-symmetric-Broyden or Bofill's formula
+//(Bofill, J. Comput. Chem. 15, 1 (1994)); TR_AUTO = exact if available, otherwise BFGS for minima and Bofill for saddle points;
+//between exact Hessians (my_hessian_recompute > 1) the updates are SR1 for minima (the exact Hessian may be indefinite) and Bofill for saddle points
+enum TR_hessian {TR_AUTO, TR_EXACT, TR_BFGS, TR_SR1, TR_PSB, TR_BOFILL};
+
+//update of a Hessian approximation B (not of its inverse) by the step s and the gradient change y; returns false if the update was skipped
+template<typename T>
+bool hessian_update(NRSMat<T> &B, const NRVec<T> &s, const NRVec<T> &y, const TR_hessian type, const T thr);
+
+template<typename T>
+class trustregion_hyperparameters : public gradient_hyperparameters<T> {
+public:
+T my_delta0; //initial trust region radius
+T my_deltamax; //maximal trust region radius
+T my_eta; //a step of a minimization is accepted if (actual reduction)/(predicted reduction) > eta
+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) {};
+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);}
+};
+
+//step of a trust-region method, used by the common driver trustregion_minimize()
+template<typename T>
+class trustregion_step {
+public:
+virtual NRVec<T> step(const NRVec<T> &g, const NRSMat<T> &B, const T delta) =0; //step within the trust region
+virtual bool saddle() const {return false;} //searching for a saddle point: steps are always taken, the radius is controlled by the quality of the model
+virtual ~trustregion_step() {};
+};
+
+//common driver of trust-region methods; achieved->my_saddle_order returns the number of negative eigenvalues of the Hessian at the result
+//(of the exact one if the problem provides it, otherwise of the approximation)
+template<typename T>
+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);
+
+
+//trust-region Newton or quasi-Newton minimization; subproblem solved exactly in the eigenbasis of B (More & Sorensen, SIAM J. Sci. Stat. Comput. 4, 553 (1983)),
+//or approximately by truncated conjugate gradients (Steihaug, SIAM J. Numer. Anal. 20, 626 (1983)), which avoids the O(n^3) diagonalization for large n
+enum TR_subproblem {TR_EIGEN, TR_STEIHAUG};
+
+template<typename T>
+class TrustRegion_hyperparameters : public trustregion_hyperparameters<T> {
+public:
+TR_subproblem my_subproblem;
+int my_maxcg; //max. no. of conjugate gradient iterations of Steihaug's method, <=0 means the dimension
+TrustRegion_hyperparameters() : my_subproblem(TR_EIGEN), my_maxcg(0) {};
+virtual void setfrom(const hyperparameters<T> &p) {const TrustRegion_hyperparameters *q=dynamic_cast<const TrustRegion_hyperparameters *>(&p); if(q) *this= *q; else trustregion_hyperparameters<T>::setfrom(p);}
+};
+
+template<typename T>
+class TrustRegion_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 ~TrustRegion_Optimizer() {};
+};
+
+
+//partitioned rational function optimization (Banerjee, Adams, Simons & Shepard, J. Phys. Chem. 89, 52 (1985); Baker, J. Comput. Chem. 7, 385 (1986)):
+//search for a stationary point with my_saddle_order negative Hessian eigenvalues, maximizing along the eigenvectors of the lowest ones,
+//minimizing along the others; saddle_order 0 gives RFO minimization
+template<typename T>
+class PRFO_hyperparameters : public trustregion_hyperparameters<T> {
+public:
+PRFO_hyperparameters() {this->my_delta0=0.3;};
+virtual void setfrom(const hyperparameters<T> &p) {const PRFO_hyperparameters *q=dynamic_cast<const PRFO_hyperparameters *>(&p); if(q) *this= *q; else trustregion_hyperparameters<T>::setfrom(p);}
+};
+
+template<typename T>
+class PRFO_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 ~PRFO_Optimizer() {};
+};
+
+
+//Levenberg-Marquardt method for least-squares problems (OptResiduals): step from (J^T J + mu D) h = -J^T r, D = I or diag(J^T J) (Marquardt scaling),
+//damping mu updated by the gain ratio (Madsen, Nielsen & Tingleff, Methods for Non-Linear Least Squares Problems, DTU 2004, alg. 3.16);
+//Jacobian from the problem or by forward differences (counted as function evaluations); convergence by |J^T r| < g_thr, by f and x changes,
+//or when the step becomes negligible at machine precision
+template<typename T>
+class LM_hyperparameters : public gradient_hyperparameters<T> {
+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) {};
+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);}
+};
+
+template<typename T>
+class LM_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 ~LM_Optimizer() {};
+};
+
+
 //Nelder-Mead downhill simplex (Lagarias et al., SIAM J. Optim. 9, 112 (1998)), derivative-free
 //convergence: f_thr is used for the spread of function values in the simplex, x_thr for its diameter, both have to be satisfied
 template<typename T>
diff --git a/t.cc b/t.cc
index 108619a..fb30c6d 100644
--- a/t.cc
+++ b/t.cc
@@ -4904,6 +4904,40 @@ const NRVec<double> &xc=x; //x is shared with fproblem.xfull, read-only access
 for(int i=0; i<n; ++i) if(frozen[i] && xopt[i]!=xc[i]) cout <<"ERROR: frozen parameter "<<i<<" changed\n";
 }
 
+if(1) //testbed of the search for saddle points: P-RFO on the Mueller-Brown potential
+{
+testproblem_mullerbrown<double> problem;
+NRVec<double> x(2);
+x[0]= -0.7; x[1]=0.6; //near the saddle point (-0.822,0.624)
+PRFO_hyperparameters<double> param;
+param.my_saddle_order=1;
+param.my_g_thr=1e-8;
+//param.my_hessian_recompute=1000; //exact Hessian only at the start, then Bofill updates
+cout << "initial " <<x<<endl;
+PRFO_Optimizer<double> o;
+hyperparameters<double> achieved;
+optresult<double> res=o.minimize(&problem,x,&param,&achieved);
+cout << "stationary point " <<x<<endl;
+cout << "converged "<<res.converged<<" V= "<<res.f<<" function evaluations "<<achieved.my_maxfunc<<" Hessian index "<<achieved.my_saddle_order<<endl;
+}
+
+if(1) //testbed of the least-squares fitting: Levenberg-Marquardt on the NIST problem Misra1a
+{
+testproblem_misra1a<double> problem;
+NRVec<double> b(2);
+for(int i=0; i<2; ++i) b[i]=problem.start(0,i);
+LM_hyperparameters<double> param;
+param.my_g_thr=1e-9;
+//param.my_scaling=true; //Marquardt scaling
+cout << "initial " <<b<<endl;
+LM_Optimizer<double> o;
+hyperparameters<double> achieved;
+optresult<double> res=o.minimize(&problem,b,&param,&achieved);
+cout << "fitted parameters " <<b<<endl;
+cout << "certified         "<<problem.certified(0)<<" "<<problem.certified(1)<<endl;
+cout << "converged "<<res.converged<<" residual sum of squares "<<2*res.f<<" certified "<<problem.rss()<<" function evaluations "<<achieved.my_maxfunc<<endl;
+}
+
 if(1) //testbed of the VQE methods, the problem provides capabilities (OptMetric, OptTrigonometric, OptHessian)
 {
 int seed;
diff --git a/testoptimize.cc b/testoptimize.cc
index dcded37..5e14d60 100644
--- a/testoptimize.cc
+++ b/testoptimize.cc
@@ -693,6 +693,124 @@ check("COBYLA with frozen parameters satisfies the constraint and keeps them",co
 }
 
 
+static void test_trustregion()
+{
+cout <<"--- trust-region methods and saddle points\n";
+{
+//secant condition B+ s = y of the Hessian updates
+srand(3);
+const int n=6;
+bool ok=true;
+double worst=0;
+const TR_hessian upd[]={TR_BFGS,TR_SR1,TR_PSB,TR_BOFILL};
+for(int t=0; t<4; ++t)
+	{
+	NRSMat<double> B(n); B.randomize(0.3); B+=2.;
+	NRVec<double> s(n), r(n); s.randomize(1.); r.randomize(0.3);
+	NRVec<double> y=B*s; y+=r;
+	if(!hessian_update(B,s,y,upd[t],1e-8)) ok=false;
+	NRVec<double> Bs=B*s; Bs-=y;
+	worst=max(worst,Bs.norm());
+	}
+check("secant condition of the BFGS, SR1, PSB and Bofill updates",ok && worst<1e-12,worst);
+}
+testproblem_mullerbrown<double> mb;
+{
+NRVec<double> x(2); x[0]= -0.3; x[1]=0.9;
+NRSMat<double> h(2), hfd(2);
+mb.hessian(x,h);
+for(int i=0; i<2; ++i) {NRVec<double> a(x), b(x), ga(2), gb(2); a.copyonwrite(); b.copyonwrite(); a[i]+=1e-6; b[i]-=1e-6; mb.funcgrad(a,&ga); mb.funcgrad(b,&gb); for(int j=0; j<=i; ++j) hfd(i,j)=(ga[j]-gb[j])/2e-6;}
+NRSMat<double> d=h-hfd;
+check("Mueller-Brown gradient and Hessian vs finite differences",fd_gradient_error(mb,x)<1e-7 && d.norm()<1e-6,d.norm());
+}
+//stationary points of the Mueller-Brown potential (Mueller & Brown 1979)
+struct {const char *name; double x0, y0, x, y, v; int k;} sp[]={
+	{"minimum (-0.558,1.442)",-0.4,1.2,-0.5582,1.4417,-146.700,0},
+	{"minimum (0.623,0.028)",0.5,0.1,0.6235,0.0280,-108.167,0},
+	{"minimum (-0.050,0.467)",0.,0.4,-0.0500,0.4667,-80.768,0},
+	{"saddle point (-0.822,0.624)",-0.7,0.6,-0.8220,0.6243,-40.665,1},
+	{"saddle point (0.212,0.293)",0.15,0.35,0.2125,0.2930,-72.249,1}};
+for(auto &c: sp) for(int m=0; m<3; ++m)
+	{
+	if(m==0 && c.k) continue;
+	NRVec<double> x(2); x[0]=c.x0; x[1]=c.y0;
+	hyperparameters<double> a;
+	bool conv;
+	double f;
+	string s;
+	if(m==0) {TrustRegion_hyperparameters<double> h; h.my_g_thr=1e-9; TrustRegion_Optimizer<double> o; f=run(o,&mb,x,h,&a,&conv); s="TR-Newton";}
+	else {PRFO_hyperparameters<double> h; h.my_g_thr=1e-9; h.my_saddle_order=c.k; if(m==2) h.my_hessian_recompute=1000; PRFO_Optimizer<double> o; f=run(o,&mb,x,h,&a,&conv); s= m==1 ? "P-RFO" : "P-RFO with Hessian updates";}
+	s+=string(" finds the ")+c.name+" with the right index";
+	check(s.c_str(),conv && abs(x[0]-c.x)<1e-4 && abs(x[1]-c.y)<1e-4 && abs(f-c.v)<1e-3 && a.my_saddle_order==c.k,f-c.v);
+	}
+{
+testproblem_vqe<double> vqe;
+TrustRegion_hyperparameters<double> h;
+TrustRegion_Optimizer<double> o;
+NRVec<double> x=startpoint(12,9);
+hyperparameters<double> a;
+bool conv;
+double f=run(o,&vqe,x,h,&a,&conv);
+check("TR-Newton with the exact Hessian on testproblem_vqe n=12",conv && f< -5.4 && a.my_maxfunc<20 && a.my_saddle_order==0,a.my_maxfunc);
+}
+{
+rosenbrock rb;
+for(int m=0; m<2; ++m)
+	{
+	TrustRegion_hyperparameters<double> h; h.my_maxstep=10000; h.my_subproblem= m ? TR_STEIHAUG : TR_EIGEN;
+	TrustRegion_Optimizer<double> o;
+	NRVec<double> x=rosenbrockstart(20);
+	bool conv;
+	double f=run(o,&rb,x,h,NULL,&conv);
+	check(m ? "TR with BFGS updates and Steihaug-CG on Rosenbrock n=20" : "TR with BFGS updates on Rosenbrock n=20",conv && f<1e-10,f);
+	}
+}
+}
+
+
+//fit of a NIST StRD problem from both starting points by Levenberg-Marquardt, certified residual sum of squares to 9 digits
+template<class P>
+static void nistfit(const char *name, const int nb)
+{
+for(int k=0; k<2; ++k)
+	{
+	P p;
+	NRVec<double> b(nb);
+	for(int i=0; i<nb; ++i) b[i]=P::start(k,i);
+	LM_hyperparameters<double> h; h.my_g_thr=1e-9; h.my_maxstep=100000;
+	LM_Optimizer<double> o;
+	bool conv;
+	double f=run(o,&p,b,h,NULL,&conv);
+	double perr=0;
+	for(int i=0; i<nb; ++i) perr=max(perr,abs(b[i]-P::certified(i))/abs(P::certified(i)));
+	string s=string("Levenberg-Marquardt reproduces the certified fit of NIST ")+name+(k ? " from start 2" : " from start 1");
+	check(s.c_str(),conv && abs(2*f-P::rss())<1e-9*P::rss() && perr<1e-5,2*f-P::rss());
+	}
+}
+
+static void test_leastsquares()
+{
+cout <<"--- least squares (Levenberg-Marquardt)\n";
+nistfit<testproblem_misra1a<double> >("Misra1a",2);
+nistfit<testproblem_mgh09<double> >("MGH09",4);
+nistfit<testproblem_thurber<double> >("Thurber (finite-difference Jacobian)",7);
+testproblem_mgh09<double> p;
+NRVec<double> b(4);
+for(int i=0; i<4; ++i) b[i]=p.start(1,i);
+NRMat<double> J(11,4), Jfd(11,4);
+NRVec<double> r(11);
+p.jacobian(b,J);
+p.residuals(b,r);
+p.jacobian_fd(b,r,Jfd);
+NRMat<double> d=J-Jfd;
+check("analytic and forward-difference Jacobians of MGH09",d.norm()<1e-6,d.norm());
+BFGS_hyperparameters<double> h; h.my_g_thr=1e-12;
+BFGS_Optimizer<double> o;
+double f=run(o,&p,b,h);
+check("other methods work on OptLeastSquares problems (BFGS on MGH09)",abs(2*f-p.rss())<1e-9*p.rss(),2*f-p.rss());
+}
+
+
 static void test_maxfunc_and_parameters()
 {
 cout <<"--- evaluation budget and hyperparameters\n";
@@ -702,9 +820,27 @@ NRVec<double> x=startpoint(6,2), xv=startpoint(8,9);
 hyperparameters<double> h; h.my_f_thr=1e-10; h.my_x_thr=1e-6; h.my_maxstep=100000;
 BFGS_Optimizer<double> bfgs; LBFGS_Optimizer<double> lbfgs; CG_Optimizer<double> cg; GDIIS_Optimizer<double> gdiis;
 NelderMead_Optimizer<double> nm; Powell_Optimizer<double> powell; Adam_Optimizer<double> adam; Rotosolve_Optimizer<double> roto; Cobyla_Optimizer<double> cobyla;
+TrustRegion_Optimizer<double> trn; PRFO_Optimizer<double> prfo; LM_Optimizer<double> lm;
+//LM evaluates the residuals directly (OptResiduals), so they are counted in the problem itself
+class countingthurber : public testproblem_thurber<double> {public: mutable int count; countingthurber() : count(0) {}; void residuals(const NRVec<double> &b, NRVec<double> &r) const {++count; testproblem_thurber<double>::residuals(b,r);}};
+{
+int bad=0;
+for(int mf=1; mf<=120; ++mf)
+	{
+	countingthurber p;
+	NRVec<double> b(7);
+	for(int i=0; i<7; ++i) b[i]=p.start(1,i);
+	h.my_maxfunc=mf;
+	hyperparameters<double> a;
+	run(lm,&p,b,h,&a);
+	if(p.count>mf || a.my_maxfunc!=p.count) ++bad;
+	}
+check("Levenberg-Marquardt (finite-difference Jacobian) respects maxfunc and reports the evaluations",bad==0,bad);
+}
 struct {const char *name; Optimizer<double> *o; const OptProblem<double> *p; NRVec<double> *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}};
+	{"Powell",&powell,&tp,&x,1},{"Adam",&adam,&tp,&x,1},{"Rotosolve",&roto,&vc,&xv,1},{"COBYLA",&cobyla,&vc,&xv,9},
+	{"TR-Newton",&trn,&vc,&xv,1},{"P-RFO",&prfo,&vc,&xv,1}};
 for(auto &c: cases)
 	{
 	int bad=maxfunc_violations(*c.o,c.p,*c.x,h,c.minfunc,120);
@@ -808,6 +944,8 @@ test_gradient_methods();
 test_derivative_free();
 test_firstorder();
 test_frozen();
+test_trustregion();
+test_leastsquares();
 test_maxfunc_and_parameters();
 test_diis();
 test_copyonwrite();
diff --git a/testoptimize.h b/testoptimize.h
index f7f1c8e..53bd17e 100644
--- a/testoptimize.h
+++ b/testoptimize.h
@@ -174,5 +174,133 @@ for(int i=0; i<nq; ++i)
 }
 
 
+
+//Mueller-Brown potential (Mueller & Brown, Theor. Chim. Acta 53, 75 (1979)), a standard 2D test surface with three minima
+//(-0.558,1.442), (0.623,0.028), (-0.050,0.467) and two saddle points (-0.822,0.624), (0.212,0.293)
+//V = sum_k A_k exp(a_k (x-x0_k)^2 + b_k (x-x0_k)(y-y0_k) + c_k (y-y0_k)^2)
+template<typename T>
+class testproblem_mullerbrown : public OptProblem<T>, public OptHessian<T> {
+public:
+static void terms(const NRVec<T> &x, const int k, T &v, T &qx, T &qy, T &a2, T &b, T &c2)
+	{
+	static const T A[]={-200,-100,-170,15}, a[]={-1,-1,-6.5,0.7}, bb[]={0,0,11,0.6}, c[]={-10,-10,-6.5,0.7}, x0[]={1,0,-0.5,-1}, y0[]={0,0.5,1.5,1};
+	const T dx=x[0]-x0[k], dy=x[1]-y0[k];
+	v=A[k]*std::exp(a[k]*dx*dx+bb[k]*dx*dy+c[k]*dy*dy);
+	qx=2*a[k]*dx+bb[k]*dy; qy=bb[k]*dx+2*c[k]*dy;
+	a2=2*a[k]; b=bb[k]; c2=2*c[k];
+	}
+virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const
+	{
+	T f=0;
+	if(grad) {grad->copyonwrite(); *grad=0.;}
+	for(int k=0; k<4; ++k)
+		{
+		T v,qx,qy,a2,b,c2;
+		terms(x,k,v,qx,qy,a2,b,c2);
+		f+=v;
+		if(grad) {(*grad)[0]+=v*qx; (*grad)[1]+=v*qy;}
+		}
+	return f;
+	}
+virtual void hessian(const NRVec<T> &x, NRSMat<T> &h) const
+	{
+	h=(T)0;
+	for(int k=0; k<4; ++k)
+		{
+		T v,qx,qy,a2,b,c2;
+		terms(x,k,v,qx,qy,a2,b,c2);
+		h(0,0)+=v*(qx*qx+a2); h(0,1)+=v*(qx*qy+b); h(1,1)+=v*(qy*qy+c2);
+		}
+	}
+virtual ~testproblem_mullerbrown() {};
+};
+
+
+
+//NIST StRD nonlinear regression test problems with certified results (https://www.itl.nist.gov/div898/strd/nls/nls_main.shtml), residuals model - data
+//start(k,i) = starting value i of "Start 1" (k=0) and "Start 2" (k=1), certified(i) = certified parameter values, rss = certified residual sum of squares
+
+//Misra1a: y = b1*(1-exp(-b2*x)), with analytic Jacobian
+template<typename T>
+class testproblem_misra1a : public OptLeastSquares<T> {
+public:
+static const T *data(const int col) //col 0: y, col 1: x
+	{
+	static const T d[2][14]={{10.07,14.73,17.94,23.93,29.61,35.18,40.02,44.82,50.76,55.05,61.01,66.40,75.47,81.78},
+				{77.6,114.9,141.1,190.8,239.9,289.0,332.8,378.4,434.8,477.3,536.8,593.1,689.1,760.0}};
+	return d[col];
+	}
+static T start(const int k, const int i) {static const T s[2][2]={{500,1e-4},{250,5e-4}}; return s[k][i];}
+static T certified(const int i) {static const T c[2]={2.3894212918E+02,5.5015643181E-04}; return c[i];}
+static T rss() {return 1.2455138894E-01;}
+virtual int nresiduals() const {return 14;}
+virtual void residuals(const NRVec<T> &b, NRVec<T> &r) const
+	{
+	r.copyonwrite();
+	for(int i=0; i<14; ++i) r[i]=b[0]*(1-std::exp(-b[1]*data(1)[i]))-data(0)[i];
+	}
+virtual bool jacobian(const NRVec<T> &b, NRMat<T> &J) const
+	{
+	J.copyonwrite();
+	for(int i=0; i<14; ++i) {const T x=data(1)[i], ex=std::exp(-b[1]*x); J(i,0)=1-ex; J(i,1)=b[0]*x*ex;}
+	return true;
+	}
+};
+
+//MGH09 (Kowalik & Osborne): y = b1*(x^2+x*b2)/(x^2+x*b3+b4), with analytic Jacobian; "Start 1" is far from the solution
+template<typename T>
+class testproblem_mgh09 : public OptLeastSquares<T> {
+public:
+static const T *data(const int col)
+	{
+	static const T d[2][11]={{1.957000E-01,1.947000E-01,1.735000E-01,1.600000E-01,8.440000E-02,6.270000E-02,4.560000E-02,3.420000E-02,3.230000E-02,2.350000E-02,2.460000E-02},
+				{4.000000E+00,2.000000E+00,1.000000E+00,5.000000E-01,2.500000E-01,1.670000E-01,1.250000E-01,1.000000E-01,8.330000E-02,7.140000E-02,6.250000E-02}};
+	return d[col];
+	}
+static T start(const int k, const int i) {static const T s[2][4]={{25,39,41.5,39},{0.25,0.39,0.415,0.39}}; return s[k][i];}
+static T certified(const int i) {static const T c[4]={1.9280693458E-01,1.9128232873E-01,1.2305650693E-01,1.3606233068E-01}; return c[i];}
+static T rss() {return 3.0750560385E-04;}
+virtual int nresiduals() const {return 11;}
+virtual void residuals(const NRVec<T> &b, NRVec<T> &r) const
+	{
+	r.copyonwrite();
+	for(int i=0; i<11; ++i) {const T x=data(1)[i]; r[i]=b[0]*(x*x+x*b[1])/(x*x+x*b[2]+b[3])-data(0)[i];}
+	}
+virtual bool jacobian(const NRVec<T> &b, NRMat<T> &J) const
+	{
+	J.copyonwrite();
+	for(int i=0; i<11; ++i)
+		{
+		const T x=data(1)[i], num=x*x+x*b[1], den=x*x+x*b[2]+b[3];
+		J(i,0)=num/den; J(i,1)=b[0]*x/den; J(i,2)= -b[0]*num*x/(den*den); J(i,3)= -b[0]*num/(den*den);
+		}
+	return true;
+	}
+};
+
+//Thurber: y = (b1 + b2*x + b3*x^2 + b4*x^3)/(1 + b5*x + b6*x^2 + b7*x^3), without Jacobian (finite differences)
+template<typename T>
+class testproblem_thurber : public OptLeastSquares<T> {
+public:
+static const T *data(const int col)
+	{
+	static const T d[2][37]={{80.574,84.248,87.264,87.195,89.076,89.608,89.868,90.101,92.405,95.854,100.696,101.060,401.672,390.724,567.534,635.316,733.054,759.087,894.206,
+				990.785,1090.109,1080.914,1122.643,1178.351,1260.531,1273.514,1288.339,1327.543,1353.863,1414.509,1425.208,1421.384,1442.962,1464.350,1468.705,1447.894,1457.628},
+				{-3.067,-2.981,-2.921,-2.912,-2.840,-2.797,-2.702,-2.699,-2.633,-2.481,-2.363,-2.322,-1.501,-1.460,-1.274,-1.212,-1.100,-1.046,-0.915,
+				-0.714,-0.566,-0.545,-0.400,-0.309,-0.109,-0.103,0.010,0.119,0.377,0.790,0.963,1.006,1.115,1.572,1.841,2.047,2.200}};
+	return d[col];
+	}
+static T start(const int k, const int i) {static const T s[2][7]={{1000,1000,400,40,0.7,0.3,0.03},{1300,1500,500,75,1,0.4,0.05}}; return s[k][i];}
+static T certified(const int i) {static const T c[7]={1.2881396800E+03,1.4910792535E+03,5.8323836877E+02,7.5416644291E+01,9.6629502864E-01,3.9797285797E-01,4.9727297349E-02}; return c[i];}
+static T rss() {return 5.6427082397E+03;}
+virtual int nresiduals() const {return 37;}
+virtual void residuals(const NRVec<T> &b, NRVec<T> &r) const
+	{
+	r.copyonwrite();
+	for(int i=0; i<37; ++i) {const T x=data(1)[i]; r[i]=(b[0]+x*(b[1]+x*(b[2]+x*b[3])))/(1+x*(b[4]+x*(b[5]+x*b[6])))-data(0)[i];}
+	}
+};
+
+
 }//namespace
 #endif /* _TESTOPTIMIZE_H_ */