diff --git a/optimize.cc b/optimize.cc
index 1424404..3d63279 100644
--- a/optimize.cc
+++ b/optimize.cc
@@ -22,7 +22,7 @@ namespace LA {
 
 
 template<typename T>
-T testproblem<T>::funcgrad(const NRVec<T> x, NRVec<T> *grad) const
+T testproblem<T>::funcgrad(const NRVec<T> &x, NRVec<T> *grad) const
 {
 int n=x.size();
 T r= n*(n-1)/2;
@@ -38,26 +38,232 @@ if(grad)
 	for(int k=0; k<n; ++k)
 		{
 		(*grad)[k] = (k+1)*x[k]*(2+0.4*x[k]*x[k]);
-		for(int i=0; i<k; ++i) (*grad)[k] += sin(x[i]*x[k]);
+		for(int i=0; i<n; ++i) if(i!=k) (*grad)[k] += x[i]*sin(x[i]*x[k]);
 		}
 	}
 return r;
 }
 
 
+//safeguarded cubic interpolation of the minimizer of phi between a and b, given function values and derivatives at both ends
 template<typename T>
-void BFGS_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par)
+static T cubic_interpolate(const T a, const T fa, const T da, const T b, const T fb, const T db)
 {
-//AI: you may suggest improvements to the general design of the class hierarchy for optimizers ( optimize.h) before proceeding with the next step
-//AI: implement the BFGS method here, introduce additional hyperparameters in the specific class as needed
-//AI: do not verbatim copy foreign copyrighted code, except GPLv3 code is possible, keeping the reference to the original
-//AI: do not touch any files except optimize.h optimize.cc and t.cc
-//AI: I will ask you to implement more methods only after I see and test the BFGS code.
+T lo= a<b ? a : b;
+T hi= a<b ? b : a;
+T safe=0.1*(hi-lo); //keep the trial point away from the interval ends
+T d1 = da + db - 3*(fa-fb)/(a-b);
+T disc = d1*d1 - da*db;
+T r;
+if(disc<0) r=(a+b)/2;
+else
+	{
+	T d2 = std::sqrt(disc); if(b<a) d2 = -d2;
+	T den = db - da + 2*d2;
+	if(den==0) r=(a+b)/2;
+	else r = b - (b-a)*(db + d2 - d1)/den;
+	}
+if(!(r>=lo+safe && r<=hi-safe)) r=(a+b)/2; //also catches NaN
+return r;
+}
+
+
+template<typename T>
+T linesearch_wolfe(const OptProblem<T> *o, NRVec<T> &x, T &f, NRVec<T> &g, const NRVec<T> &p, T alpha1, const T c1, const T c2, const T alphamax, const int maxiter, int &nfunc, const int verbose)
+{
+const int n=x.size();
+const NRVec<T> x0(x);
+const T f0=f;
+const T d0= g*p;
+if(d0>=0) {if(verbose) std::cout <<"linesearch: not a descent direction\n"; return 0;}
+if(alpha1>alphamax) alpha1=alphamax;
+
+//data of trial points; vectors are always freshly allocated, so that assignments (which share data) are safe
+T alpha0=0, f_alpha0=f0, d_alpha0=d0; NRVec<T> x_alpha0(x0), g_alpha0(g); //previous trial point
+T alo=0,flo=0,dlo=0,ahi=0,fhi=0,dhi=0; NRVec<T> xlo,glo; //bracketing interval for the zoom phase, lo is the best point
+bool bracketed=false;
+int it;
+
+//bracketing phase
+for(it=0; it<maxiter; ++it)
+	{
+	NRVec<T> xt = x0 + p*alpha1;
+	NRVec<T> gt(n);
+	T ft=o->funcgrad(xt,&gt); ++nfunc;
+	T dt= gt*p;
+	if(verbose>1) std::cout <<"linesearch bracket "<<it<<" alpha= "<<alpha1<<" f= "<<ft<<" dphi= "<<dt<<std::endl;
+	if(ft > f0 + c1*alpha1*d0 || (it>0 && ft>=f_alpha0))
+		{
+		alo=alpha0; flo=f_alpha0; dlo=d_alpha0; xlo=x_alpha0; glo=g_alpha0;
+		ahi=alpha1; fhi=ft; dhi=dt;
+		bracketed=true;
+		break;
+		}
+	if(std::abs(dt) <= -c2*d0 || alpha1>=alphamax) {x=xt; f=ft; g=gt; return alpha1;} //strong Wolfe satisfied or cannot go further
+	if(dt>=0)
+		{
+		alo=alpha1; flo=ft; dlo=dt; xlo=xt; glo=gt;
+		ahi=alpha0; fhi=f_alpha0; dhi=d_alpha0;
+		bracketed=true;
+		break;
+		}
+	alpha0=alpha1; f_alpha0=ft; d_alpha0=dt; x_alpha0=xt; g_alpha0=gt;
+	alpha1 *= 2; if(alpha1>alphamax) alpha1=alphamax;
+	}
+
+//zoom phase, invariant: alo satisfies sufficient decrease and has the lowest f found so far, dlo*(ahi-alo)<0
+if(bracketed) for(++it; it<maxiter; ++it)
+	{
+	T a=cubic_interpolate(alo,flo,dlo,ahi,fhi,dhi);
+	NRVec<T> xt = x0 + p*a;
+	NRVec<T> gt(n);
+	T ft=o->funcgrad(xt,&gt); ++nfunc;
+	T dt= gt*p;
+	if(verbose>1) std::cout <<"linesearch zoom "<<it<<" alpha= "<<a<<" f= "<<ft<<" dphi= "<<dt<<std::endl;
+	if(ft > f0 + c1*a*d0 || ft>=flo)
+		{
+		ahi=a; fhi=ft; dhi=dt;
+		}
+	else
+		{
+		if(std::abs(dt) <= -c2*d0) {x=xt; f=ft; g=gt; return a;}
+		if(dt*(ahi-alo)>=0) {ahi=alo; fhi=flo; dhi=dlo;}
+		alo=a; flo=ft; dlo=dt; xlo=xt; glo=gt;
+		}
+	if(std::abs(ahi-alo) <= std::numeric_limits<T>::epsilon()*std::abs(alo)) break; //interval collapsed
+	}
+
+//failure - accept the best point satisfying sufficient decrease, if any
+if(!bracketed && alpha0>0) {alo=alpha0; flo=f_alpha0; xlo=x_alpha0; glo=g_alpha0;} //bracketing ran out of evaluations
+if(alo>0)
+	{
+	if(verbose) std::cout <<"linesearch: strong Wolfe conditions not satisfied, accepting sufficient decrease step "<<alo<<std::endl;
+	x=xlo; f=flo; g=glo; return alo;
+	}
+if(verbose) std::cout <<"linesearch failed\n";
+return 0;
+}
+
+
+//BFGS update of the inverse Hessian approximation H, cf. Nocedal & Wright eq. 6.17
+//H+ = (I - rho s y^T) H (I - rho y s^T) + rho s s^T, rho=1/(y^T s)
+//expanded: H+ = H + (rho^2 y^T H y + rho) s s^T - rho (s (Hy)^T + Hy s^T) = H + s u^T + u s^T with u = (rho^2 y^T H y + rho)/2 s - rho Hy
+template<typename T>
+static void BFGS_update(NRSMat<T> &H, const NRVec<T> &s, const NRVec<T> &y, const T sy)
+{
+T rho = 1/sy;
+NRVec<T> Hy = H*y;
+T yHy = y*Hy;
+NRVec<T> u = s*((rho*rho*yHy + rho)/2);
+u.axpy(-rho,Hy);
+H.rank2update((T)1,s,u);
+}
+
+
+template<typename T>
+optresult<T> BFGS_Optimizer<T>::minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par, hyperparameters<T> *achieved)
+{
+//collect parameters, those not provided by the caller are taken from defaults
+BFGS_hyperparameters<T> hp;
+if(par)
+	{
+	static_cast<hyperparameters<T> &>(hp) = *par;
+	const gradient_hyperparameters<T> *gpar = dynamic_cast<const gradient_hyperparameters<T> *>(par);
+	if(gpar) static_cast<gradient_hyperparameters<T> &>(hp) = *gpar;
+	const BFGS_hyperparameters<T> *bpar = dynamic_cast<const BFGS_hyperparameters<T> *>(par);
+	if(bpar) hp = *bpar;
+	}
+const int verbose=hp.my_verbose;
+
+const int n=x.size();
+NRVec<T> g(n);
+T f=o->funcgrad(x,&g);
+int nfunc=1;
+NRSMat<T> H(n);
+H=1.;
+bool Hinit=true; //H is the initial guess, rescale it after the first step
+bool converged=false;
+T last_df=0, last_dx=0;
+
+int iter;
+for(iter=0; iter<hp.my_maxstep; ++iter)
+	{
+	T gnorm=g.norm();
+	if(verbose) std::cout <<"BFGS iter "<<iter<<" f= "<<f<<" |g|= "<<gnorm<<std::endl;
+	if(gnorm<hp.my_g_thr) {converged=true; break;}
+	if(hp.my_maxfunc>0 && nfunc>=hp.my_maxfunc) {if(verbose) std::cout <<"BFGS: maximum number of function evaluations reached\n"; break;}
+
+	NRVec<T> p = H*g;
+	p *= (T)-1;
+	if(p*g >= 0) //H lost positive definiteness numerically, restart
+		{
+		if(verbose) std::cout <<"BFGS: not a descent direction, resetting inverse Hessian\n";
+		H=1.; Hinit=true;
+		p = g*(T)-1;
+		}
+
+	NRVec<T> xold(x), gold(g);
+	T fold=f;
+	T pnorm=p.norm();
+	T alpha1 = Hinit ? 1/pnorm : 1; //on the first step (or after restart) try unit step length in x, otherwise the Newton step
+	T alphamax = hp.my_maxsteplength/pnorm;
+	int maxls=hp.my_maxlinesearch;
+	if(hp.my_maxfunc>0 && hp.my_maxfunc-nfunc<maxls) maxls=hp.my_maxfunc-nfunc;
+	T alpha=linesearch_wolfe(o,x,f,g,p,alpha1,hp.my_c1,hp.my_c2,alphamax,maxls,nfunc,verbose);
+	if(alpha==0)
+		{
+		if(hp.my_maxfunc>0 && nfunc>=hp.my_maxfunc) continue; //reported at the top of the loop
+		if(Hinit) {if(verbose) std::cout <<"BFGS: line search failed along steepest descent, giving up\n"; break;}
+		if(verbose) std::cout <<"BFGS: line search failed, resetting inverse Hessian\n";
+		H=1.; Hinit=true;
+		continue;
+		}
+
+	NRVec<T> s = x-xold;
+	NRVec<T> y = g-gold;
+	T sy = s*y;
+	T snorm=s.norm();
+	T ynorm=y.norm();
+	if(verbose>1) std::cout <<"BFGS step alpha= "<<alpha<<" |s|= "<<snorm<<" s.y= "<<sy<<std::endl;
+	if(sy > hp.my_sy_thr*snorm*ynorm)
+		{
+		if(Hinit) {H = (T)(sy/(ynorm*ynorm)); Hinit=false;} //scale initial guess, Nocedal & Wright eq. 6.20
+		BFGS_update(H,s,y,sy);
+		}
+	else if(verbose) std::cout <<"BFGS: curvature condition violated, skipping update\n";
+
+	last_df=std::abs(f-fold); last_dx=snorm;
+	if(last_df<hp.my_f_thr && last_dx<hp.my_x_thr) {converged=true; ++iter; break;}
+	}
+
+if(verbose)
+	{
+	if(converged) std::cout <<"BFGS converged in "<<iter<<" iterations, "<<nfunc<<" function evaluations, f= "<<f<<" |g|= "<<g.norm()<<std::endl;
+	else std::cout <<"BFGS NOT converged after "<<iter<<" iterations, "<<nfunc<<" function evaluations, f= "<<f<<" |g|= "<<g.norm()<<std::endl;
+	}
+
+if(achieved)
+	{
+	BFGS_hyperparameters<T> ach(hp);
+	ach.my_maxstep=iter;
+	ach.my_maxfunc=nfunc;
+	ach.my_f_thr=last_df;
+	ach.my_x_thr=last_dx;
+	ach.my_g_thr=g.norm();
+	*achieved = static_cast<hyperparameters<T> &>(ach);
+	gradient_hyperparameters<T> *gach = dynamic_cast<gradient_hyperparameters<T> *>(achieved);
+	if(gach) *gach = static_cast<gradient_hyperparameters<T> &>(ach);
+	BFGS_hyperparameters<T> *bach = dynamic_cast<BFGS_hyperparameters<T> *>(achieved);
+	if(bach) *bach = ach;
+	}
+
+return optresult<T>(f,converged);
 }
 
 
 
 //force instantiation
+template double linesearch_wolfe(const OptProblem<double> *o, NRVec<double> &x, double &f, NRVec<double> &g, const NRVec<double> &p, double alpha1, const double c1, const double c2, const double alphamax, const int maxiter, int &nfunc, const int verbose);
 template class testproblem<double>;
 template class BFGS_Optimizer<double>;
 
diff --git a/optimize.h b/optimize.h
index d35c10f..64439ab 100644
--- a/optimize.h
+++ b/optimize.h
@@ -37,7 +37,7 @@ namespace LA {
 template<typename T>
 class OptProblem {
 public:
-virtual T funcgrad(const NRVec<T> x, NRVec<T> *grad) const =0; //evaluate function and optionally also gradient at point x
+virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const =0; //evaluate function and optionally also gradient at point x (if grad!=NULL)
 virtual ~OptProblem() {};
 };
 
@@ -45,7 +45,7 @@ virtual ~OptProblem() {};
 template<typename T>
 class testproblem : public OptProblem<T> {
 public:
-virtual T funcgrad(const NRVec<T> x, NRVec<T> *grad) const;
+virtual T funcgrad(const NRVec<T> &x, NRVec<T> *grad) const;
 virtual ~testproblem() {};
 };
 
@@ -53,43 +53,77 @@ virtual ~testproblem() {};
 //base classes for optimization methods
 
 //this is a collection of parameters, which can be augmented by method-specific additions, to be used as a base class
+//methods obtain their specific parameters via dynamic_cast, missing ones are taken from the default constructor
 template<typename T>
 class hyperparameters {
 public:
 int my_maxstep; //maximum number of iterations
+int my_maxfunc; //maximum number of function evaluations, <=0 means unlimited
 T my_f_thr; //threshhold on function change in a single step
 T my_x_thr; //threshold on x norm change in a single step
-
-virtual int maxstep() const {return my_maxstep;}
-virtual T f_thr() const {return my_f_thr;}
-virtual T x_thr() const {return my_x_thr;}
-virtual T g_thr() const =0;
+int my_verbose; //0=silent, 1=per-iteration info, 2=more detailed info
+hyperparameters() : my_maxstep(200), my_maxfunc(0), my_f_thr(0), my_x_thr(0), my_verbose(1) {};
 virtual ~hyperparameters(){};
 };
 
+//parameters common to methods using the gradient
+template<typename T>
+class gradient_hyperparameters : public hyperparameters<T> {
+public:
+T my_g_thr; //gradient norm convergence threshold
+gradient_hyperparameters() : my_g_thr(1e-6) {};
+};
+
+
+//result of the minimization, converts implicitly to the function value
+template<typename T>
+class optresult {
+public:
+T f; //final function value
+bool converged;
+optresult(const T ff, const bool c) : f(ff), converged(c) {};
+operator T() const {return f;}
+};
+
 
 //base class for optimization methods, virtual functions to allow a pointer to derived class
 template<typename T>
 class Optimizer {
 public:
-virtual void minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par) = 0;
+//par==NULL means default parameters
+//if achieved!=NULL, the values actually reached are returned there: my_maxstep = no. of iterations done, my_maxfunc = no. of function evaluations done, thresholds = last values of the convergence criteria
+//(|f change| and |x change| in the last step, 0 if no step was done; |gradient| for gradient methods), the other parameters as they were used
+virtual optresult<T> minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par=NULL, hyperparameters<T> *achieved=NULL) = 0;
 virtual ~Optimizer() {}; 
 };
 
 
+//auxiliary routines
+
+//line search along direction p satisfying strong Wolfe conditions (cf. Nocedal & Wright, Numerical Optimization, 2nd ed., alg. 3.5 and 3.6)
+//on input x,f,g at the starting point, on successful output x,f,g at the new point; returns step length or 0 on failure (x,f,g then unchanged)
+//maxiter = max. number of function evaluations, nfunc is incremented by the number of evaluations done
+template<typename T>
+T linesearch_wolfe(const OptProblem<T> *o, NRVec<T> &x, T &f, NRVec<T> &g, const NRVec<T> &p, T alpha1, const T c1, const T c2, const T alphamax, const int maxiter, int &nfunc, const int verbose);
+
+
 //specific methods
 
 template<typename T>
-class BFGS_hyperparameters : public hyperparameters<T> {
+class BFGS_hyperparameters : public gradient_hyperparameters<T> {
 public:
-T my_g_thr;
-virtual T g_thr() const {return my_g_thr;} //gradient magnitude convergence threshold
+T my_c1; //sufficient decrease (Armijo) parameter of Wolfe conditions, 0<c1<c2<1
+T my_c2; //curvature parameter of Wolfe conditions
+int my_maxlinesearch; //max. number of function evaluations in one line search
+T my_maxsteplength; //max. norm of a single step in x
+T my_sy_thr; //skip the update if s.y <= sy_thr*|s|*|y| (curvature condition violated numerically)
+BFGS_hyperparameters() : my_c1(1e-4), my_c2(0.9), my_maxlinesearch(30), my_maxsteplength(1e10), my_sy_thr(1e-12) {};
 };
 
 template<typename T>
 class BFGS_Optimizer : public Optimizer<T> {
 public:
-virtual void minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par);
+virtual optresult<T> minimize(const OptProblem<T> *o, NRVec<T> &x, const hyperparameters<T> *par=NULL, hyperparameters<T> *achieved=NULL);
 virtual ~BFGS_Optimizer() {};
 };
 
diff --git a/t.cc b/t.cc
index de9d2fd..252877f 100644
--- a/t.cc
+++ b/t.cc
@@ -4850,8 +4850,10 @@ param.my_x_thr=1e-3;
 param.my_g_thr=1e-3;
 cout << "initial " <<x<<endl;
 BFGS_Optimizer<double> o;
-o.minimize(&problem,x,&param);
+BFGS_hyperparameters<double> achieved;
+optresult<double> res=o.minimize(&problem,x,&param,&achieved);
 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;
 double r=problem.funcgrad(x,NULL);
 cout <<r<<endl;
 }