diff --git a/fastlib/trunk/contrib/tqlong/optimization/CMakeLists.txt b/fastlib/trunk/contrib/tqlong/optimization/CMakeLists.txt
index bd2692330c..65864068f6 100644
--- a/fastlib/trunk/contrib/tqlong/optimization/CMakeLists.txt
+++ b/fastlib/trunk/contrib/tqlong/optimization/CMakeLists.txt
@@ -5,6 +5,7 @@ add_executable(optim_tqlong_test
    test1.cpp
    FunctionTemplate
    NelderMead
+   GradientDescent
 )
 # link dependencies of executable
 target_link_libraries(optim_tqlong_test
diff --git a/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.cpp b/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.cpp
index 1973845af9..8323521280 100644
--- a/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.cpp
+++ b/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.cpp
@@ -4,6 +4,7 @@
 
 double optim::LengthEuclidianSquare::CalculateValue(const Vector &x)
 {
+  DEBUG_ASSERT(x.length() == dim);
   double s = 0;
   for (int i = 0; i < x.length(); i++)
     s += x[i]*x[i];
@@ -12,13 +13,14 @@ double optim::LengthEuclidianSquare::CalculateValue(const Vector &x)
 
 void optim::LengthEuclidianSquare::CalculateGradient(const Vector &x, Vector &gradient)
 {
+  DEBUG_ASSERT(x.length() == dim && gradient.length() == dim);
   DEBUG_ASSERT(x.length() == gradient.length());
   la::ScaleOverwrite(2.0, x, &gradient);
 }
 
 void optim::LengthEuclidianSquare::CalculateHessian(const Vector &x, Matrix &hessian)
 {
-  DEBUG_ASSERT(x.length() == hessian.n_rows() && x.length() == hessian.n_cols());
+  DEBUG_ASSERT(x.length() == hessian.n_rows() && x.length() == hessian.n_cols() && x.length() == dim);
   hessian.SetAll(0.0);
   for (int i = 0; i < hessian.n_rows(); i++)
     hessian.ref(i, i) = 2.0;
diff --git a/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.h b/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.h
index 96b1d0b0b4..3f208869f3 100644
--- a/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.h
+++ b/fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.h
@@ -14,16 +14,26 @@ BEGIN_OPTIM_NAMESPACE;
 
 /*************************************************************
   * A function template, implement CalculateXXXXX methods
-  * 0-order smooth function: CalculateValue
+  * Function value         : CalculateValue
   * 1-order smooth function: CalculateValue, CalculateGradient
   * 2-order smooth function: CalculateValue, CalculateGradient, CalculateHessian
+  * Required:
+      int dimension()
+      void Init(variable_type*)
+      variable_type (typedef)
   *************************************************************/
 
 class LengthEuclidianSquare {
+  int dim;
 public:
+  typedef Vector variable_type;                   // required
+  int dimension() { return dim; }                 // required
+  void Init(Vector* x) { x->Init(dim); }          // required
   double CalculateValue(const Vector& x);
   void CalculateGradient(const Vector& x, Vector& gradient);
   void CalculateHessian(const Vector& x, Matrix& hessian);
+public:
+  LengthEuclidianSquare(int dim_ = 2) : dim(dim_) {}
 };
 
 END_OPTIM_NAMESPACE;
diff --git a/fastlib/trunk/contrib/tqlong/optimization/GradientDescent.h b/fastlib/trunk/contrib/tqlong/optimization/GradientDescent.h
new file mode 100644
index 0000000000..88d133b5ae
--- /dev/null
+++ b/fastlib/trunk/contrib/tqlong/optimization/GradientDescent.h
@@ -0,0 +1,234 @@
+#ifndef GRADIENTDESCENT_H
+#define GRADIENTDESCENT_H
+
+#include <fastlib/fastlib.h>
+
+#ifndef BEGIN_OPTIM_NAMESPACE
+#define BEGIN_OPTIM_NAMESPACE namespace optim {
+#endif
+#ifndef END_OPTIM_NAMESPACE
+#define END_OPTIM_NAMESPACE }
+#endif
+
+BEGIN_OPTIM_NAMESPACE;
+
+/**********************************************************************
+  Implement 1-order Gradient Descent optimization method with Wolfe line search method
+  Function: CalculateValue(const variable_type&)
+            CalculateGradient(const variable_type&, variable_type& gradient)
+  Function::variable_type:
+    implement la::AddExpert(double, const variable_type&, variable_type*)
+              la::ScaleOverwrite(double, const variable_type&, variable_type*);
+              la::LengthEuclidean(const variable_type&)
+              la::Dot(const variable_type&, const variable_type&)
+              variable_type.CopyValues(const variable_type&)
+  Parameter:
+    General params
+      maxIter, rTol, aTol
+    Specific for Gradient Descent & Wolfe line search
+      c1, c2, beta
+**********************************************************************/
+template<typename Function> class GradientDescent
+{
+public:
+  typedef Function function_type;
+  typedef typename Function::variable_type variable_type;
+  typedef double* OptimizationParameters;
+    //maximum # of iterations  : maxIter = (int) param[0];
+    //relative tolerance       : rTol = param[1];
+    //absolute tolerance       : aTol = param[2];
+    //c1 : Wolfe 1st condition : c1 = param[3];
+    //c2 : Wolfe 2nd condition : c2 = param[4];
+    //beta : scale parameter   : beta = param[5];
+  struct HistoryRecord {
+    int iter;
+    int n_evals;
+    int n_grads;
+    double best_val;
+    double residual;
+    OT_DEF(HistoryRecord) {
+      OT_MY_OBJECT(iter);
+      OT_MY_OBJECT(n_evals);
+      OT_MY_OBJECT(n_grads);
+      OT_MY_OBJECT(best_val);
+      OT_MY_OBJECT(residual);
+    }
+  public:
+    HistoryRecord(int iter_, int n_evals_, int n_grads_, double best_val_, double residual_)
+      : iter(iter_), n_evals(n_evals_), n_grads(n_grads_), best_val(best_val_), residual(residual_) {}
+  };
+protected:
+  static double default_gradient_descent_parameter[6]; // = {100, 0.001, 0.01, 1e-4, 0.9, 0.5};
+public:
+  GradientDescent(function_type& f_, OptimizationParameters param_ = default_gradient_descent_parameter);
+  void setParam(OptimizationParameters param);
+  void setX0(const variable_type& x0_);
+  double optimize(variable_type& sol);
+  void printHistory();
+
+  ArrayList<HistoryRecord> history;
+protected:
+  function_type& f;
+  variable_type x0;
+
+  // parameters
+  OptimizationParameters param;
+  int maxIter;
+  double aTol, rTol;
+  double c1, c2, beta;
+
+  // progress
+  int iter;
+  int n_evals;
+  int n_grads;
+  double best_val;
+  double residual;
+
+  double WolfeStep(const variable_type& x, double val_x, const variable_type& grad, const variable_type& p);
+  double CalculateValue(const variable_type& x);
+  void CalculateGradient(const variable_type& x, variable_type& grad);
+  void recordProgress();
+  void printProgress();
+};
+
+template<typename F>
+double GradientDescent<F>::CalculateValue(const variable_type& x) {
+  double val = f.CalculateValue(x);
+  n_evals++;
+  return val;
+}
+
+template<typename F>
+void GradientDescent<F>::CalculateGradient(const variable_type& x, variable_type& grad) {
+  f.CalculateGradient(x, grad);
+  n_grads++;
+}
+
+template<typename F>
+void GradientDescent<F>::recordProgress() {
+  history.PushBackCopy(HistoryRecord(iter, n_evals, n_grads, best_val, residual));
+  printProgress();
+}
+
+template<typename F>
+void GradientDescent<F>::printProgress() {
+  if (history.size() == 0) return;
+  printf("iter = %d n_evals = %d n_grads = %d best_val = %f residual = %f\n",
+         history.back().iter, history.back().n_evals, history.back().n_grads,
+         history.back().best_val, history.back().residual);
+}
+
+template<typename F>
+void GradientDescent<F>::printHistory() {
+  printf("History:\n");
+  for (int i = 0; i < history.size(); i++) {
+    printf("iter = %d n_evals = %d n_grads = %d best_val = %f residual = %f\n",
+           history[i].iter, history[i].n_evals, history[i].n_grads, history[i].best_val, history[i].residual);
+  }
+}
+
+template<typename F>
+double GradientDescent<F>::default_gradient_descent_parameter[6] = {100, 0.001, 0.01, 1e-4, 0.9, 0.5};
+
+template<typename F>
+GradientDescent<F>::GradientDescent(function_type &f_, OptimizationParameters param_)
+  : f(f_), param(param_) {
+  setParam(param);
+  f.Init(&x0);
+  history.Init();
+}
+
+template<typename F>
+void GradientDescent<F>::setParam(OptimizationParameters param) {
+  maxIter = (int) param[0];
+  rTol = param[1];
+  aTol = param[2];
+  c1 = param[3];
+  c2 = param[4];
+  beta = param[5];
+}
+
+template<typename F>
+void GradientDescent<F>::setX0(const variable_type& x0_) {
+  x0.CopyValues(x0_);
+}
+
+template<typename F>
+double GradientDescent<F>::optimize(variable_type &sol) {
+  history.Clear();
+  iter = 0;
+  n_evals = 0;
+  n_grads = 0;
+
+  variable_type x; // the current search variable
+  variable_type grad, d; // the current search direction
+
+  f.Init(&x);      // initialize search variable
+  f.Init(&d);      // and direction
+  f.Init(&grad);      // and direction
+
+  x.CopyValues(x0);
+  double val = CalculateValue(x);
+  CalculateGradient(x, grad); // Calculate gradient
+  double r0 = la::LengthEuclidean(grad);
+
+  sol.CopyValues(x0);
+  best_val = val;
+  residual = r0;
+  recordProgress();
+  for (iter = 1; iter < maxIter; iter++) {
+    // Calculate search direction: negative gradient
+    la::ScaleOverwrite(-1.0, grad, &d);
+    // Calculate step size by Wolfe's conditions
+    double lambda = WolfeStep(x, val, grad, d);
+    // Calculate new variable
+    la::AddExpert(lambda, d, &x);
+    // Update best value
+    val = CalculateValue(x);
+    if (val < best_val) {
+      best_val = val;
+      sol.CopyValues(x);
+    }
+    // Calculate gradient and check termination condition
+    CalculateGradient(x, grad);
+    residual = la::LengthEuclidean(grad);
+    recordProgress();
+    if (residual < rTol*r0+aTol) break;
+  }
+  return best_val;
+}
+
+template<typename F>
+double GradientDescent<F>::WolfeStep(const variable_type& x, double val_x, const variable_type& grad, const variable_type& p) {
+  double lambda = 1.0/beta, val_xp;
+  double dot_grad_p;
+  variable_type x_p, grad_p;
+  f.Init(&x_p);
+  f.Init(&grad_p);
+
+  dot_grad_p = la::Dot(grad, p);
+  while (1) {
+    lambda *= beta;
+    if (lambda < 1e-10) {
+      printf("Line search results in a too small step size, try increasing c2 (param[4]).\n");
+      break;
+    }
+    x_p.CopyValues(x);
+    la::AddExpert(lambda, p, &x_p);     // x_p = x + lambda*p
+
+    val_xp = CalculateValue(x_p);     // f(x_p)
+
+    // first Wolfe condition
+    if (val_xp - val_x <= c1*lambda*dot_grad_p) {  // f(x_p) - f(x) <= c1 * lambda * <grad, p>
+      // second Wolfe condition
+      CalculateGradient(x_p, grad_p);   // grad f(x_p)
+      double dot_grad_x_p = la::Dot(grad_p, p);
+      if (dot_grad_x_p >= c2*dot_grad_p) break;          // <p,grad f(x_p)> >= c2 * <p, grad>
+    }
+  }
+  return lambda;
+}
+
+END_OPTIM_NAMESPACE;
+
+#endif // GRADIENTDESCENT_H
diff --git a/fastlib/trunk/contrib/tqlong/optimization/NelderMead.h b/fastlib/trunk/contrib/tqlong/optimization/NelderMead.h
index c654470ee0..e1f6e53df9 100644
--- a/fastlib/trunk/contrib/tqlong/optimization/NelderMead.h
+++ b/fastlib/trunk/contrib/tqlong/optimization/NelderMead.h
@@ -24,20 +24,26 @@ struct OptimizationParameter {
 BEGIN_OPTIM_NAMESPACE;
 
 /**********************************************************************
-  Implement 0-order Nelder-Mead optimization method
-  Variable: implement la::AddExpert(alpha, const Variable& X, Variable* Y) and Copy() for this type
-  Function: implement CalculateValue(const Variable&)
+  Implement 0-order Nelder-Mead simplex optimization method
+  Function: implement CalculateValue(const variable_type&)
+  Function::variable_type:
+    implement la::AddExpert(double, const variable_type&, variable_type*)
+              Copy(const variable_type&), CopyValues(const variable_type&)
+              for this type
   Parameter:
-    alpha : Reflection  (1.0)
-    gamma : Expansion   (2.0)
-    rho   : Contraction (0.5)
-    sigma : Reduction   (0.5)
+    General params
+      maxIter, rTol, aTol
+    Specific for Nelder-Mead
+      alpha : Reflection  (1.0)
+      gamma : Expansion   (2.0)
+      rho   : Contraction (0.5)
+      sigma : Reduction   (0.5)
 **********************************************************************/
-template<typename Variable, typename Function> class NelderMead
+template<typename Function> class NelderMead
 {
 public:
   typedef Function function_type;
-  typedef Variable variable_type;
+  typedef typename Function::variable_type variable_type;
 public:
   NelderMead(function_type& f_, OptimizationParameter param_ = OptimizationParameter(100, 0.01, 0.01));
   void setParam(double alpha_, double gamma_, double rho_, double sigma_) {
@@ -64,33 +70,33 @@ protected:
   void updateCenter(const variable_type& x1, const variable_type& x2);
 };
 
-template<typename V, typename F>
-NelderMead<V, F>::NelderMead(function_type& f_, OptimizationParameter param_)
+template<typename F>
+NelderMead<F>::NelderMead(function_type& f_, OptimizationParameter param_)
   : f(f_), param(param_)
 {
   memory.Init();
   val.Init();
+  f.Init(&center);
   alpha = 1.5;
   gamma = 2.0;
   rho = 0.5;
   sigma = 0.5;
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::add(const variable_type &x)
+template<typename F>
+void NelderMead<F>::add(const variable_type &x)
 {
   double v = f.CalculateValue(x);
   add(x, v);
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::addSeed(const ArrayList<variable_type>& vX) {
-  for (int i = 0; i < vX.size(); i++)
-    add(vX[i]);
+template<typename F>
+void NelderMead<F>::addSeed(const ArrayList<variable_type>& vX) {
+  for (int i = 0; i < vX.size(); i++) add(vX[i]);
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::add(const variable_type &x, double v)
+template<typename F>
+void NelderMead<F>::add(const variable_type &x, double v)
 {
   double pos = findPos(v);
   if (pos < memory.size()) {
@@ -104,19 +110,19 @@ void NelderMead<V, F>::add(const variable_type &x, double v)
   updateCenter(x);
 }
 
-template<typename V, typename F>
-int NelderMead<V, F>::findPos(double v)
+template<typename F>
+int NelderMead<F>::findPos(double v)
 {
   int i = 0;
   while (i < val.size() && val[i] < v) i++;
   return i;
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::updateCenter(const variable_type &x)
+template<typename F>
+void NelderMead<F>::updateCenter(const variable_type &x)
 {
   if (memory.size() <= 1)
-    center.Copy(x);
+    center.CopyValues(x);
   else {
     variable_type tmp;
     tmp.Copy(center);
@@ -126,8 +132,8 @@ void NelderMead<V, F>::updateCenter(const variable_type &x)
   }
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::updateCenter(const variable_type &x1, const variable_type &x2)
+template<typename F>
+void NelderMead<F>::updateCenter(const variable_type &x1, const variable_type &x2)
 {
   int n = memory.size();
   DEBUG_ASSERT(n > 0);
@@ -135,8 +141,8 @@ void NelderMead<V, F>::updateCenter(const variable_type &x1, const variable_type
   la::AddExpert(1.0/n, x2, &center);
 }
 
-template<typename V, typename F>
-double NelderMead<V, F>::optimize(variable_type &sol)
+template<typename F>
+double NelderMead<F>::optimize(variable_type &sol)
 {
   int n = memory.size();
   DEBUG_ASSERT(n >= 3);
@@ -195,8 +201,8 @@ double NelderMead<V, F>::optimize(variable_type &sol)
   return best_val;
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::replace(int pos, const variable_type &x, double v)
+template<typename F>
+void NelderMead<F>::replace(int pos, const variable_type &x, double v)
 {
   int n = memory.size();
   DEBUG_ASSERT(pos < n);
@@ -218,8 +224,8 @@ void NelderMead<V, F>::replace(int pos, const variable_type &x, double v)
   val[pos] = v;
 }
 
-template<typename V, typename F>
-void NelderMead<V, F>::updateAll()
+template<typename F>
+void NelderMead<F>::updateAll()
 {
   int n = memory.size();
   DEBUG_ASSERT(n > 0);
diff --git a/fastlib/trunk/contrib/tqlong/optimization/test1.cpp b/fastlib/trunk/contrib/tqlong/optimization/test1.cpp
index c76d9bb869..df22440c38 100644
--- a/fastlib/trunk/contrib/tqlong/optimization/test1.cpp
+++ b/fastlib/trunk/contrib/tqlong/optimization/test1.cpp
@@ -3,6 +3,7 @@
 
 #include "FunctionTemplate.h"
 #include "NelderMead.h"
+#include "GradientDescent.h"
 
 using namespace std;
 using namespace optim;
@@ -27,20 +28,22 @@ void testNelderMead(fx_module* module) {
   cout << "Nelder-Mead test ..." << endl;
 
   int n_seed = fx_param_int(module, "seed", 10);
+  LengthEuclidianSquare f(2);
+  NelderMead<optim::LengthEuclidianSquare> algo(f);
+
+  // Seeding
   ArrayList<Vector> vX;
   vX.Init();
   for (int i = 0; i < n_seed; i++) {
-    Vector x; x.Init(2);
+    Vector x;
+    f.Init(&x);
     x[0] = 10+10*(double)rand()/RAND_MAX;
     x[1] = 20+10*(double)rand()/RAND_MAX;
     vX.PushBackCopy(x);
   }
-  LengthEuclidianSquare f;
-
-  NelderMead<Vector, optim::LengthEuclidianSquare> algo(f);
-
   algo.addSeed(vX);
 
+  // Optimization
   Vector sol;
   sol.Init(2);
   double v = algo.optimize(sol);
@@ -51,11 +54,38 @@ void testNelderMead(fx_module* module) {
   cout << "Nelder-Mead test succeeded." << endl;
 }
 
+void testGradientDescent(fx_module* module) {
+  cout << "GradientDescent test ..." << endl;
+
+  double param[] = {100, 0.00001, 0.001, 1e-4, 0.9, 0.4};
+  LengthEuclidianSquare f(2);
+  GradientDescent<optim::LengthEuclidianSquare> algo(f, param);
+
+  // Seeding
+  Vector x0;
+  f.Init(&x0);
+  x0[0] = 10+10*(double)rand()/RAND_MAX;
+  x0[1] = 20+10*(double)rand()/RAND_MAX;
+
+  // Optimization
+  Vector sol;
+  f.Init(&sol);
+  algo.setX0(x0);
+  double v = algo.optimize(sol);
+
+  cout << "Best value = " << v << endl;
+  ot::Print(sol, "Solution", stdout);
+  //algo.printHistory();
+
+  cout << "GradientDescent test succeeded." << endl;
+}
+
 int main(int argc, char** argv) {
   fx_module* root = fx_init(argc, argv, &optimization_doc);
   cout << "Optimization tests" << endl;
 
-  testNelderMead(root);
+  //testNelderMead(root);
+  testGradientDescent(root);
 
   fx_done(root);
   return 0;