added scg

minimize.cc

@@ -330,8 +330,7 @@ void minimize(ColumnVector& x, const EvalFunction& func){
 
     sort(v.begin(),v.end(),pair_comparer()); //double bracks constructs a temporary object
 
-//      cout << "how="<< how << endl;
-//      cout << "en "<<v[0].first<<endl;
+
   } //closing while
  
   x=v[0].second;
@@ -339,24 +338,112 @@ void minimize(ColumnVector& x, const EvalFunction& func){
   
 }
 
-//  int main(int argc,char *argv[])
-//  {
-
-//    ColumnVector x(4);
-//    x <<3 <<1<<14<<9;;
-//    Matrix hess4,hess2,hess1;
 
-//    minimize(x);
-//    cerr<<"params"<<endl;
-//    cerr <<x<<endl;
+void scg(ColumnVector& x,const gEvalFunction& func){
+  int niters=100; //default maybe make an options option
+  
+  int fevals=0;
+  int gevals=0;
+  int nparams=x.Nrows();
+  float sigma0 = 1.0e-4;
+  float  fold=func.evaluate(x); fevals++;
+  float fnow=0,fnew=0;
+  ColumnVector gradold=func.g_evaluate(x);gevals++;
+  ColumnVector gradnew=gradold;
+  ColumnVector d=-gradnew;       // search direction
+  ColumnVector xplus,xnew,gplus;
+  bool success=true;
+  int nsuccess=0;
+  float lambda=1.0;
+  float   lambdamin = 1.0e-15; 
+  float lambdamax = 1.0e100;			
+  int j = 1;	
+  float mu=0,kappa=0,sigma=0,gamma=0,alpha=0,delta=0,Delta,beta=0;
+  float eps=1e-16;
+  // main loop..
+  while(j<niters){
+
+    if(success){
+      mu=(d.t()*gradnew).AsScalar();
+      if(mu >= 0){
+	d=-gradnew;
+	mu=(d.t()*gradnew).AsScalar();
+      }
+      kappa=(d.t()*d).AsScalar();
+      if(kappa<eps){
+	break;
+      }
+      sigma=sigma0/sqrt(kappa);
+      xplus = x + sigma*d;
+      gplus=func.g_evaluate(xplus);gevals++;
+      gamma = (d.t()*(gplus - gradnew)).AsScalar()/sigma;
+    }
+    delta = gamma + lambda*kappa;
+    if (delta <= 0){ 
+      delta = lambda*kappa;
+      lambda = lambda - gamma/kappa;
+    }
+    alpha = - mu/delta;
+    
 
-//    hess4=hessian(x,0.0001,4);
-//    hess2=hessian(x,0.0001,2);
-//    hess1=hessian(x,0.0001,1);
-//    cerr<<"Hessian"<<endl;
-//    cerr<<"order1"<<endl<<hess1<<endl<<"order2"<<endl<<hess2<<endl<<endl<<"order4"<<endl<<hess4;
+     xnew = x + alpha*d;
+     fnew=func.evaluate(xnew);fevals++;
+     Delta = 2*(fnew - fold)/(alpha*mu);
+     if (Delta  >= 0){
+       success = true;
+       nsuccess = nsuccess + 1;
+       x = xnew;
+       fnow = fnew;}
+     else{
+       success = false;
+       fnow = fold;
+     }
+     
+     if (success == 1){
+     
+
+       //Test for termination...
+       if(0==1){
+	 break;
+       }
+       else{
+	 fold = fnew;
+	 gradold = gradnew;
+	 gradnew=func.g_evaluate(x);gevals++;
+	 if ((gradnew.t()*gradnew).AsScalar() == 0){
+	   break;
+	 }
+       }
+       
+     }
+     
+     if (Delta < 0.25){
+       //   lambda = min(4.0*lambda, lambdamax);
+       lambda=4.0*lambda<lambdamax ? 4.0*lambda : lambdamax;
+     }
+     if (Delta > 0.75){
+       //lambda = max(0.5*lambda, lambdamin);
+       lambda = 0.5*lambda > lambdamin ? 0.5*lambda : lambdamin;
+     }
+     
+
+     if (nsuccess == nparams){
+       d = -gradnew;
+       nsuccess = 0;
+     }
+     else{
+       if (success == 1){	 
+	 beta = ((gradold - gradnew).t()*gradnew).AsScalar()/mu;
+	 d = beta*d - gradnew;
+       }
+     }
+     j++; 
+  }
   
-//  }
+}
+
+
+
 
 }
 

minimize.h

@@ -42,11 +42,10 @@ public:
 };
 
 class EvalFunction
-{
+{//Function where gradient is not analytic (or you are too lazy to work it out)
 public:
-  
   EvalFunction(){}
-  virtual float evaluate(const ColumnVector& x) const = 0;
+  virtual float evaluate(const ColumnVector& x) const = 0; //evaluate the function
   virtual ~EvalFunction(){};
 
 private:
@@ -55,17 +54,37 @@ private:
   EvalFunction(const EvalFunction&);
 };
 
-float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);
+class gEvalFunction : public EvalFunction
+{//Function where gradient is analytic (required for scg)
+public:
+  gEvalFunction() : EvalFunction(){}
+  // evaluate is inherited from EvalFunction
+  // virtual float evaluate(const ColumnVector& x) const = 0; //evaluate the function  
+  virtual ColumnVector g_evaluate(const ColumnVector& x) const = 0; //evaluate the gradient
+  virtual ~gEvalFunction(){};
+
+private:
+    
+  const gEvalFunction& operator=(gEvalFunction& par);
+  gEvalFunction(const gEvalFunction&);
+};
+
+
+float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff derivative
 
-float diff2(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);
+float diff2(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff 2nd derivative
 
-float diff2(const ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);
+float diff2(const ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);// finite diff cross derivative
+
+ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff hessian
 
-ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);
 
 void minimize(ColumnVector& x, const EvalFunction& func);
+void scg(ColumnVector& x, const gEvalFunction& func);
+}
+   
+
 
-}   
 
 #endif