git-svn-id: file:///home/svn/opt_utilities@129 ed2142bd-67ad-457f-ba7c-d818d4011675

4 files changed, 1963 insertions, 0 deletions

diff --git a/methods/lbfgs/arithmetic_ansi.h b/methods/lbfgs/arithmetic_ansi.h
new file mode 100644
index 0000000..c58c98a
--- /dev/null
+++ b/methods/lbfgs/arithmetic_ansi.h
@@ -0,0 +1,133 @@
+/*
+ *      ANSI C implementation of vector operations.
+ *
+ * Copyright (c) 2007-2010 Naoaki Okazaki
+ * All rights reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+/* $Id: arithmetic_ansi.h 65 2010-01-29 12:19:16Z naoaki $ */
+
+#include <stdlib.h>
+#include <memory.h>
+
+#if     LBFGS_FLOAT == 32 && LBFGS_IEEE_FLOAT
+#define fsigndiff(x, y) (((*(uint32_t*)(x)) ^ (*(uint32_t*)(y))) & 0x80000000U)
+#else
+#define fsigndiff(x, y) (*(x) * (*(y) / fabs(*(y))) < 0.)
+#endif/*LBFGS_IEEE_FLOAT*/
+
+inline static void* vecalloc(size_t size)
+{
+    void *memblock = malloc(size);
+    if (memblock) {
+        memset(memblock, 0, size);
+    }
+    return memblock;
+}
+
+inline static void vecfree(void *memblock)
+{
+    free(memblock);
+}
+
+inline static void vecset(lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
+{
+    int i;
+    
+    for (i = 0;i < n;++i) {
+        x[i] = c;
+    }
+}
+
+inline static void veccpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
+{
+    int i;
+
+    for (i = 0;i < n;++i) {
+        y[i] = x[i];
+    }
+}
+
+inline static void vecncpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
+{
+    int i;
+
+    for (i = 0;i < n;++i) {
+        y[i] = -x[i];
+    }
+}
+
+inline static void vecadd(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
+{
+    int i;
+
+    for (i = 0;i < n;++i) {
+        y[i] += c * x[i];
+    }
+}
+
+inline static void vecdiff(lbfgsfloatval_t *z, const lbfgsfloatval_t *x, const lbfgsfloatval_t *y, const int n)
+{
+    int i;
+
+    for (i = 0;i < n;++i) {
+        z[i] = x[i] - y[i];
+    }
+}
+
+inline static void vecscale(lbfgsfloatval_t *y, const lbfgsfloatval_t c, const int n)
+{
+    int i;
+
+    for (i = 0;i < n;++i) {
+        y[i] *= c;
+    }
+}
+
+inline static void vecmul(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
+{
+    int i;
+
+    for (i = 0;i < n;++i) {
+        y[i] *= x[i];
+    }
+}
+
+inline static void vecdot(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const lbfgsfloatval_t *y, const int n)
+{
+    int i;
+    *s = 0.;
+    for (i = 0;i < n;++i) {
+        *s += x[i] * y[i];
+    }
+}
+
+inline static void vec2norm(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n)
+{
+    vecdot(s, x, x, n);
+    *s = (lbfgsfloatval_t)sqrt(*s);
+}
+
+inline static void vec2norminv(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n)
+{
+    vec2norm(s, x, n);
+    *s = (lbfgsfloatval_t)(1.0 / *s);
+}
diff --git a/methods/lbfgs/lbfgs.cpp b/methods/lbfgs/lbfgs.cpp
new file mode 100644
index 0000000..d7b5494
--- /dev/null
+++ b/methods/lbfgs/lbfgs.cpp
@@ -0,0 +1,1276 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "arithmetic_ansi.h"
+
+
+//#define min2(a, b)      ((a) <= (b) ? (a) : (b))
+template<typename T>
+static T min2(const T& a,const T& b)
+{
+  return (a) <= (b) ? (a) : (b);
+}
+//#define max2(a, b)      ((a) >= (b) ? (a) : (b))
+template<typename T>
+static T max2(const T& a,const T& b)
+{
+  return (a) >= (b) ? (a) : (b);
+}
+//#define max3(a, b, c)   max2(max2((a), (b)), (c));
+template<typename T>
+static T max3(const T& a,const T& b,const T& c)
+{
+  return max2(max2((a), (b)), (c));
+}
+
+
+struct tag_callback_data {
+  int n;
+  void *instance;
+  lbfgs_evaluate_t proc_evaluate;
+  lbfgs_progress_t proc_progress;
+};
+typedef struct tag_callback_data callback_data_t;
+
+struct tag_iteration_data {
+  lbfgsfloatval_t alpha;
+  lbfgsfloatval_t *s;     /* [n] */
+  lbfgsfloatval_t *y;     /* [n] */
+  lbfgsfloatval_t ys;     /* vecdot(y, s) */
+};
+typedef struct tag_iteration_data iteration_data_t;
+
+static const lbfgs_parameter_t _defparam = {
+  6, 1e-5, 0, 1e-5,
+  0, LBFGS_LINESEARCH_DEFAULT, 40,
+  1e-20, 1e20, 1e-4, 0.9, 0.9, 1.0e-16,
+  0.0, 0, -1,
+};
+
+/* Forward function declarations. */
+
+typedef int (*line_search_proc)(
+				int n,
+				lbfgsfloatval_t *x,
+				lbfgsfloatval_t *f,
+				lbfgsfloatval_t *g,
+				lbfgsfloatval_t *s,
+				lbfgsfloatval_t *stp,
+				const lbfgsfloatval_t* xp,
+				const lbfgsfloatval_t* gp,
+				lbfgsfloatval_t *wa,
+				callback_data_t *cd,
+				const lbfgs_parameter_t *param
+				);
+
+static int line_search_backtracking(
+				    int n,
+				    lbfgsfloatval_t *x,
+				    lbfgsfloatval_t *f,
+				    lbfgsfloatval_t *g,
+				    lbfgsfloatval_t *s,
+				    lbfgsfloatval_t *stp,
+				    const lbfgsfloatval_t* xp,
+				    const lbfgsfloatval_t* gp,
+				    lbfgsfloatval_t *wa,
+				    callback_data_t *cd,
+				    const lbfgs_parameter_t *param
+				    );
+
+static int line_search_backtracking_owlqn(
+					  int n,
+					  lbfgsfloatval_t *x,
+					  lbfgsfloatval_t *f,
+					  lbfgsfloatval_t *g,
+					  lbfgsfloatval_t *s,
+					  lbfgsfloatval_t *stp,
+					  const lbfgsfloatval_t* xp,
+					  const lbfgsfloatval_t* gp,
+					  lbfgsfloatval_t *wp,
+					  callback_data_t *cd,
+					  const lbfgs_parameter_t *param
+					  );
+
+static int line_search_morethuente(
+				   int n,
+				   lbfgsfloatval_t *x,
+				   lbfgsfloatval_t *f,
+				   lbfgsfloatval_t *g,
+				   lbfgsfloatval_t *s,
+				   lbfgsfloatval_t *stp,
+				   const lbfgsfloatval_t* xp,
+				   const lbfgsfloatval_t* gp,
+				   lbfgsfloatval_t *wa,
+				   callback_data_t *cd,
+				   const lbfgs_parameter_t *param
+				   );
+
+static int update_trial_interval(
+				 lbfgsfloatval_t *x,
+				 lbfgsfloatval_t *fx,
+				 lbfgsfloatval_t *dx,
+				 lbfgsfloatval_t *y,
+				 lbfgsfloatval_t *fy,
+				 lbfgsfloatval_t *dy,
+				 lbfgsfloatval_t *t,
+				 lbfgsfloatval_t *ft,
+				 lbfgsfloatval_t *dt,
+				 const lbfgsfloatval_t tmin,
+				 const lbfgsfloatval_t tmax,
+				 int *brackt
+				 );
+
+static lbfgsfloatval_t owlqn_x1norm(
+				    const lbfgsfloatval_t* x,
+				    const int start,
+				    const int n
+				    );
+
+static void owlqn_pseudo_gradient(
+				  lbfgsfloatval_t* pg,
+				  const lbfgsfloatval_t* x,
+				  const lbfgsfloatval_t* g,
+				  const int n,
+				  const lbfgsfloatval_t c,
+				  const int start,
+				  const int end
+				  );
+
+static void owlqn_project(
+			  lbfgsfloatval_t* d,
+			  const lbfgsfloatval_t* sign,
+			  const int start,
+			  const int end
+			  );
+
+
+static lbfgsfloatval_t* lbfgs_malloc(int n)
+{
+  return (lbfgsfloatval_t*)vecalloc(sizeof(lbfgsfloatval_t) * n);
+}
+
+static void lbfgs_free(lbfgsfloatval_t *x)
+{
+  vecfree(x);
+}
+
+static void lbfgs_parameter_init(lbfgs_parameter_t *param)
+{
+  memcpy(param, &_defparam, sizeof(*param));
+}
+
+static int lbfgs(
+	  int n,
+	  lbfgsfloatval_t *x,
+	  lbfgsfloatval_t *ptr_fx,
+	  lbfgs_evaluate_t proc_evaluate,
+	  lbfgs_progress_t proc_progress,
+	  void *instance,
+	  lbfgs_parameter_t *_param
+	  )
+{
+  int ret;
+  int i, j, k, ls, end, bound;
+  lbfgsfloatval_t step;
+  
+  /* Constant parameters and their default values. */
+  lbfgs_parameter_t param = (_param != NULL) ? (*_param) : _defparam;
+  const int m = param.m;
+  
+  lbfgsfloatval_t *xp = NULL;
+  lbfgsfloatval_t *g = NULL, *gp = NULL, *pg = NULL;
+  lbfgsfloatval_t *d = NULL, *w = NULL, *pf = NULL;
+  iteration_data_t *lm = NULL, *it = NULL;
+  lbfgsfloatval_t ys, yy;
+  lbfgsfloatval_t xnorm, gnorm, beta;
+  lbfgsfloatval_t fx = 0.;
+  lbfgsfloatval_t rate = 0.;
+  line_search_proc linesearch = line_search_morethuente;
+  
+  /* Construct a callback data. */
+  callback_data_t cd;
+  cd.n = n;
+  cd.instance = instance;
+  cd.proc_evaluate = proc_evaluate;
+  cd.proc_progress = proc_progress;
+  
+  /* Check the input parameters for errors. */
+  if (n <= 0) {
+    return LBFGSERR_INVALID_N;
+  }
+  
+  if (param.epsilon < 0.) {
+    return LBFGSERR_INVALID_EPSILON;
+  }
+  if (param.past < 0) {
+    return LBFGSERR_INVALID_TESTPERIOD;
+  }
+  if (param.delta < 0.) {
+    return LBFGSERR_INVALID_DELTA;
+  }
+  if (param.min_step < 0.) {
+    return LBFGSERR_INVALID_MINSTEP;
+  }
+  if (param.max_step < param.min_step) {
+    return LBFGSERR_INVALID_MAXSTEP;
+  }
+  if (param.ftol < 0.) {
+    return LBFGSERR_INVALID_FTOL;
+  }
+  if (param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE ||
+      param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE) {
+    if (param.wolfe <= param.ftol || 1. <= param.wolfe) {
+      return LBFGSERR_INVALID_WOLFE;
+    }
+  }
+  if (param.gtol < 0.) {
+    return LBFGSERR_INVALID_GTOL;
+  }
+  if (param.xtol < 0.) {
+    return LBFGSERR_INVALID_XTOL;
+  }
+  if (param.max_linesearch <= 0) {
+    return LBFGSERR_INVALID_MAXLINESEARCH;
+  }
+  if (param.orthantwise_c < 0.) {
+    return LBFGSERR_INVALID_ORTHANTWISE;
+  }
+  if (param.orthantwise_start < 0 || n < param.orthantwise_start) {
+    return LBFGSERR_INVALID_ORTHANTWISE_START;
+  }
+  if (param.orthantwise_end < 0) {
+    param.orthantwise_end = n;
+  }
+  if (n < param.orthantwise_end) {
+    return LBFGSERR_INVALID_ORTHANTWISE_END;
+  }
+  if (param.orthantwise_c != 0.) {
+    switch (param.linesearch) {
+    case LBFGS_LINESEARCH_BACKTRACKING:
+      linesearch = line_search_backtracking_owlqn;
+      break;
+    default:
+      /* Only the backtracking method is available. */
+      return LBFGSERR_INVALID_LINESEARCH;
+    }
+  } else {
+    switch (param.linesearch) {
+    case LBFGS_LINESEARCH_MORETHUENTE:
+      linesearch = line_search_morethuente;
+      break;
+    case LBFGS_LINESEARCH_BACKTRACKING_ARMIJO:
+    case LBFGS_LINESEARCH_BACKTRACKING_WOLFE:
+    case LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE:
+      linesearch = line_search_backtracking;
+      break;
+    default:
+      return LBFGSERR_INVALID_LINESEARCH;
+    }
+  }
+  
+  /* Allocate working space. */
+  xp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+  g = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+  gp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+  d = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+  w = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+  if (xp == NULL || g == NULL || gp == NULL || d == NULL || w == NULL) {
+    ret = LBFGSERR_OUTOFMEMORY;
+    goto lbfgs_exit;
+  }
+  
+  if (param.orthantwise_c != 0.) {
+    /* Allocate working space for OW-LQN. */
+    pg = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+    if (pg == NULL) {
+      ret = LBFGSERR_OUTOFMEMORY;
+      goto lbfgs_exit;
+    }
+  }
+  
+  /* Allocate limited memory storage. */
+  lm = (iteration_data_t*)vecalloc(m * sizeof(iteration_data_t));
+  if (lm == NULL) {
+    ret = LBFGSERR_OUTOFMEMORY;
+    goto lbfgs_exit;
+  }
+  
+  /* Initialize the limited memory. */
+  for (i = 0;i < m;++i) {
+    it = &lm[i];
+    it->alpha = 0;
+    it->ys = 0;
+    it->s = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+    it->y = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+    if (it->s == NULL || it->y == NULL) {
+      ret = LBFGSERR_OUTOFMEMORY;
+      goto lbfgs_exit;
+    }
+  }
+  
+  /* Allocate an array for storing previous values of the objective function. */
+  if (0 < param.past) {
+    pf = (lbfgsfloatval_t*)vecalloc(param.past * sizeof(lbfgsfloatval_t));
+  }
+  
+  /* Evaluate the function value and its gradient. */
+  fx = cd.proc_evaluate(cd.instance, x, g, cd.n, 0);
+  if (0. != param.orthantwise_c) {
+    /* Compute the L1 norm of the variable and add it to the object value. */
+    xnorm = owlqn_x1norm(x, param.orthantwise_start, param.orthantwise_end);
+    fx += xnorm * param.orthantwise_c;
+    owlqn_pseudo_gradient(
+			  pg, x, g, n,
+			  param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
+			  );
+  }
+  
+  /* Store the initial value of the objective function. */
+  if (pf != NULL) {
+    pf[0] = fx;
+  }
+  
+  /*
+    Compute the direction;
+    we assume the initial hessian matrix H_0 as the identity matrix.
+  */
+  if (param.orthantwise_c == 0.) {
+    vecncpy(d, g, n);
+  } else {
+    vecncpy(d, pg, n);
+  }
+  
+  /*
+    Make sure that the initial variables are not a minimizer.
+  */
+  vec2norm(&xnorm, x, n);
+  if (param.orthantwise_c == 0.) {
+    vec2norm(&gnorm, g, n);
+  } else {
+    vec2norm(&gnorm, pg, n);
+  }
+  if (xnorm < 1.0) xnorm = 1.0;
+  if (gnorm / xnorm <= param.epsilon) {
+    ret = LBFGS_ALREADY_MINIMIZED;
+    goto lbfgs_exit;
+  }
+  
+  /* Compute the initial step:
+     step = 1.0 / sqrt(vecdot(d, d, n))
+  */
+  vec2norminv(&step, d, n);
+  
+  k = 1;
+  end = 0;
+  for (;;) {
+    /* Store the current position and gradient vectors. */
+    veccpy(xp, x, n);
+    veccpy(gp, g, n);
+    
+    /* Search for an optimal step. */
+    if (param.orthantwise_c == 0.) {
+      ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, &param);
+    } else {
+      ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, &param);
+      owlqn_pseudo_gradient(
+			    pg, x, g, n,
+			    param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
+			    );
+    }
+    if (ls < 0) {
+      /* Revert to the previous point. */
+      veccpy(x, xp, n);
+      veccpy(g, gp, n);
+      ret = ls;
+      goto lbfgs_exit;
+    }
+    
+    /* Compute x and g norms. */
+    vec2norm(&xnorm, x, n);
+    if (param.orthantwise_c == 0.) {
+      vec2norm(&gnorm, g, n);
+    } else {
+      vec2norm(&gnorm, pg, n);
+    }
+    
+    /* Report the progress. */
+    if (cd.proc_progress) {
+      if (ret = cd.proc_progress(cd.instance, x, g, fx, xnorm, gnorm, step, cd.n, k, ls)) {
+	goto lbfgs_exit;
+      }
+    }
+    
+    /*
+      Convergence test.
+      The criterion is given by the following formula:
+      |g(x)| / \max(1, |x|) < \epsilon
+    */
+    if (xnorm < 1.0) xnorm = 1.0;
+    if (gnorm / xnorm <= param.epsilon) {
+      /* Convergence. */
+      ret = LBFGS_SUCCESS;
+      break;
+    }
+    
+    /*
+      Test for stopping criterion.
+      The criterion is given by the following formula:
+      (f(past_x) - f(x)) / f(x) < \delta
+    */
+    if (pf != NULL) {
+      /* We don't test the stopping criterion while k < past. */
+      if (param.past <= k) {
+	/* Compute the relative improvement from the past. */
+	rate = (pf[k % param.past] - fx) / fx;
+	
+	/* The stopping criterion. */
+	if (rate < param.delta) {
+	  ret = LBFGS_STOP;
+	  break;
+	}
+      }
+      
+      /* Store the current value of the objective function. */
+      pf[k % param.past] = fx;
+    }
+    
+    if (param.max_iterations != 0 && param.max_iterations < k+1) {
+      /* Maximum number of iterations. */
+      ret = LBFGSERR_MAXIMUMITERATION;
+      break;
+    }
+    
+    /*
+      Update vectors s and y:
+      s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}.
+      y_{k+1} = g_{k+1} - g_{k}.
+    */
+    it = &lm[end];
+    vecdiff(it->s, x, xp, n);
+    vecdiff(it->y, g, gp, n);
+    
+    /*
+      Compute scalars ys and yy:
+      ys = y^t \cdot s = 1 / \rho.
+      yy = y^t \cdot y.
+      Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor).
+    */
+    vecdot(&ys, it->y, it->s, n);
+    vecdot(&yy, it->y, it->y, n);
+    it->ys = ys;
+    
+    /*
+      Recursive formula to compute dir = -(H \cdot g).
+      This is described in page 779 of:
+      Jorge Nocedal.
+      Updating Quasi-Newton Matrices with Limited Storage.
+      Mathematics of Computation, Vol. 35, No. 151,
+      pp. 773--782, 1980.
+    */
+    bound = (m <= k) ? m : k;
+    ++k;
+    end = (end + 1) % m;
+    
+    /* Compute the steepest direction. */
+    if (param.orthantwise_c == 0.) {
+      /* Compute the negative of gradients. */
+      vecncpy(d, g, n);
+    } else {
+      vecncpy(d, pg, n);
+    }
+    
+    j = end;
+    for (i = 0;i < bound;++i) {
+      j = (j + m - 1) % m;    /* if (--j == -1) j = m-1; */
+      it = &lm[j];
+      /* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */
+      vecdot(&it->alpha, it->s, d, n);
+      it->alpha /= it->ys;
+      /* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */
+      vecadd(d, it->y, -it->alpha, n);
+    }
+    
+    vecscale(d, ys / yy, n);
+    
+    for (i = 0;i < bound;++i) {
+      it = &lm[j];
+      /* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamma_{i}. */
+      vecdot(&beta, it->y, d, n);
+      beta /= it->ys;
+      /* \gamma_{i+1} = \gamma_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */
+      vecadd(d, it->s, it->alpha - beta, n);
+      j = (j + 1) % m;        /* if (++j == m) j = 0; */
+    }
+    
+    /*
+      Constrain the search direction for orthant-wise updates.
+    */
+    if (param.orthantwise_c != 0.) {
+      for (i = param.orthantwise_start;i < param.orthantwise_end;++i) {
+	if (d[i] * pg[i] >= 0) {
+	  d[i] = 0;
+	}
+      }
+    }
+    
+    /*
+      Now the search direction d is ready. We try step = 1 first.
+    */
+    step = 1.0;
+  }
+  
+ lbfgs_exit:
+  /* Return the final value of the objective function. */
+  if (ptr_fx != NULL) {
+    *ptr_fx = fx;
+  }
+  
+  vecfree(pf);
+  
+  /* Free memory blocks used by this function. */
+  if (lm != NULL) {
+    for (i = 0;i < m;++i) {
+      vecfree(lm[i].s);
+      vecfree(lm[i].y);
+    }
+    vecfree(lm);
+  }
+  vecfree(pg);
+  vecfree(w);
+  vecfree(d);
+  vecfree(gp);
+  vecfree(g);
+  vecfree(xp);
+  
+  return ret;
+}
+
+
+
+static int line_search_backtracking(
+				    int n,
+				    lbfgsfloatval_t *x,
+				    lbfgsfloatval_t *f,
+				    lbfgsfloatval_t *g,
+				    lbfgsfloatval_t *s,
+				    lbfgsfloatval_t *stp,
+				    const lbfgsfloatval_t* xp,
+				    const lbfgsfloatval_t* gp,
+				    lbfgsfloatval_t *wp,
+				    callback_data_t *cd,
+				    const lbfgs_parameter_t *param
+				    )
+{
+  int ret = 0, count = 0;
+  lbfgsfloatval_t width, dg, norm = 0.;
+  lbfgsfloatval_t finit, dginit = 0., dgtest;
+  const lbfgsfloatval_t dec = 0.5, inc = 2.1;
+  
+  /* Check the input parameters for errors. */
+  if (*stp <= 0.) {
+    return LBFGSERR_INVALIDPARAMETERS;
+  }
+  
+  /* Compute the initial gradient in the search direction. */
+  vecdot(&dginit, g, s, n);
+  
+  /* Make sure that s points to a descent direction. */
+  if (0 < dginit) {
+    return LBFGSERR_INCREASEGRADIENT;
+  }
+  
+  /* The initial value of the objective function. */
+  finit = *f;
+  dgtest = param->ftol * dginit;
+  
+  for (;;) {
+    veccpy(x, xp, n);
+    vecadd(x, s, *stp, n);
+    
+    /* Evaluate the function and gradient values. */
+    *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp);
+    
+    ++count;
+    
+    if (*f > finit + *stp * dgtest) {
+      width = dec;
+    } else {
+      /* The sufficient decrease condition (Armijo condition). */
+      if (param->linesearch == LBFGS_LINESEARCH_BACKTRACKING_ARMIJO) {
+	/* Exit with the Armijo condition. */
+	return count;
+      }
+      
+      /* Check the Wolfe condition. */
+      vecdot(&dg, g, s, n);
+      if (dg < param->wolfe * dginit) {
+	width = inc;
+      } else {
+	if(param->linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE) {
+	  /* Exit with the regular Wolfe condition. */
+	  return count;
+	}
+	
+	/* Check the strong Wolfe condition. */
+	if(dg > -param->wolfe * dginit) {
+	  width = dec;
+	} else {
+	  /* Exit with the strong Wolfe condition. */
+	  return count;
+	}
+      }
+    }
+    
+    if (*stp < param->min_step) {
+      /* The step is the minimum value. */
+      return LBFGSERR_MINIMUMSTEP;
+    }
+    if (*stp > param->max_step) {
+      /* The step is the maximum value. */
+      return LBFGSERR_MAXIMUMSTEP;
+    }
+    if (param->max_linesearch <= count) {
+      /* Maximum number of iteration. */
+      return LBFGSERR_MAXIMUMLINESEARCH;
+    }
+    
+    (*stp) *= width;
+  }
+}
+
+
+
+static int line_search_backtracking_owlqn(
+					  int n,
+					  lbfgsfloatval_t *x,
+					  lbfgsfloatval_t *f,
+					  lbfgsfloatval_t *g,
+					  lbfgsfloatval_t *s,
+					  lbfgsfloatval_t *stp,
+					  const lbfgsfloatval_t* xp,
+					  const lbfgsfloatval_t* gp,
+					  lbfgsfloatval_t *wp,
+					  callback_data_t *cd,
+					  const lbfgs_parameter_t *param
+					  )
+{
+  int i, ret = 0, count = 0;
+  lbfgsfloatval_t width = 0.5, norm = 0.;
+  lbfgsfloatval_t finit = *f, dgtest;
+  
+  /* Check the input parameters for errors. */
+  if (*stp <= 0.) {
+    return LBFGSERR_INVALIDPARAMETERS;
+  }
+  
+  /* Choose the orthant for the new point. */
+  for (i = 0;i < n;++i) {
+    wp[i] = (xp[i] == 0.) ? -gp[i] : xp[i];
+  }
+  
+  for (;;) {
+    /* Update the current point. */
+    veccpy(x, xp, n);
+    vecadd(x, s, *stp, n);
+    
+    /* The current point is projected onto the orthant. */
+    owlqn_project(x, wp, param->orthantwise_start, param->orthantwise_end);
+    
+    /* Evaluate the function and gradient values. */
+    *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp);
+    
+    /* Compute the L1 norm of the variables and add it to the object value. */
+    norm = owlqn_x1norm(x, param->orthantwise_start, param->orthantwise_end);
+    *f += norm * param->orthantwise_c;
+    
+    ++count;
+    
+    dgtest = 0.;
+    for (i = 0;i < n;++i) {
+      dgtest += (x[i] - xp[i]) * gp[i];
+    }
+    
+    if (*f <= finit + param->ftol * dgtest) {
+      /* The sufficient decrease condition. */
+      return count;
+    }
+    
+    if (*stp < param->min_step) {
+      /* The step is the minimum value. */
+      return LBFGSERR_MINIMUMSTEP;
+    }
+    if (*stp > param->max_step) {
+      /* The step is the maximum value. */
+      return LBFGSERR_MAXIMUMSTEP;
+    }
+    if (param->max_linesearch <= count) {
+      /* Maximum number of iteration. */
+      return LBFGSERR_MAXIMUMLINESEARCH;
+    }
+    
+    (*stp) *= width;
+  }
+}
+
+
+
+static int line_search_morethuente(
+				   int n,
+				   lbfgsfloatval_t *x,
+				   lbfgsfloatval_t *f,
+				   lbfgsfloatval_t *g,
+				   lbfgsfloatval_t *s,
+				   lbfgsfloatval_t *stp,
+				   const lbfgsfloatval_t* xp,
+				   const lbfgsfloatval_t* gp,
+				   lbfgsfloatval_t *wa,
+				   callback_data_t *cd,
+				   const lbfgs_parameter_t *param
+				   )
+{
+  int count = 0;
+  int brackt, stage1, uinfo = 0;
+  lbfgsfloatval_t dg;
+  lbfgsfloatval_t stx, fx, dgx;
+  lbfgsfloatval_t sty, fy, dgy;
+  lbfgsfloatval_t fxm, dgxm, fym, dgym, fm, dgm;
+  lbfgsfloatval_t finit, ftest1, dginit, dgtest;
+  lbfgsfloatval_t width, prev_width;
+  lbfgsfloatval_t stmin, stmax;
+  
+  /* Check the input parameters for errors. */
+  if (*stp <= 0.) {
+    return LBFGSERR_INVALIDPARAMETERS;
+  }
+  
+  /* Compute the initial gradient in the search direction. */
+  vecdot(&dginit, g, s, n);
+  
+  /* Make sure that s points to a descent direction. */
+  if (0 < dginit) {
+    return LBFGSERR_INCREASEGRADIENT;
+  }
+  
+  /* Initialize local variables. */
+  brackt = 0;
+  stage1 = 1;
+  finit = *f;
+  dgtest = param->ftol * dginit;
+  width = param->max_step - param->min_step;
+  prev_width = 2.0 * width;
+  
+  /*
+    The variables stx, fx, dgx contain the values of the step,
+    function, and directional derivative at the best step.
+    The variables sty, fy, dgy contain the value of the step,
+    function, and derivative at the other endpoint of
+    the interval of uncertainty.
+    The variables stp, f, dg contain the values of the step,
+    function, and derivative at the current step.
+  */
+  stx = sty = 0.;
+  fx = fy = finit;
+  dgx = dgy = dginit;
+  
+  for (;;) {
+    /*
+      Set the minimum and maximum steps to correspond to the
+      present interval of uncertainty.
+    */
+    if (brackt) {
+      stmin = min2(stx, sty);
+      stmax = max2(stx, sty);
+    } else {
+      stmin = stx;
+      stmax = *stp + 4.0 * (*stp - stx);
+    }
+    
+    /* Clip the step in the range of [stpmin, stpmax]. */
+    if (*stp < param->min_step) *stp = param->min_step;
+    if (param->max_step < *stp) *stp = param->max_step;
+    
+    /*
+      If an unusual termination is to occur then let
+      stp be the lowest point obtained so far.
+    */
+    if ((brackt && ((*stp <= stmin || stmax <= *stp) || param->max_linesearch <= count + 1 || uinfo != 0)) || (brackt && (stmax - stmin <= param->xtol * stmax))) {
+      *stp = stx;
+    }
+	
+    /*
+      Compute the current value of x:
+      x <- x + (*stp) * s.
+    */
+    veccpy(x, xp, n);
+    vecadd(x, s, *stp, n);
+	
+    /* Evaluate the function and gradient values. */
+    *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp);
+    vecdot(&dg, g, s, n);
+	
+    ftest1 = finit + *stp * dgtest;
+    ++count;
+	
+    /* Test for errors and convergence. */
+    if (brackt && ((*stp <= stmin || stmax <= *stp) || uinfo != 0)) {
+      /* Rounding errors prevent further progress. */
+      return LBFGSERR_ROUNDING_ERROR;
+    }
+    if (*stp == param->max_step && *f <= ftest1 && dg <= dgtest) {
+      /* The step is the maximum value. */
+      return LBFGSERR_MAXIMUMSTEP;
+    }
+    if (*stp == param->min_step && (ftest1 < *f || dgtest <= dg)) {
+      /* The step is the minimum value. */
+      return LBFGSERR_MINIMUMSTEP;
+    }
+    if (brackt && (stmax - stmin) <= param->xtol * stmax) {
+      /* Relative width of the interval of uncertainty is at most xtol. */
+      return LBFGSERR_WIDTHTOOSMALL;
+    }
+    if (param->max_linesearch <= count) {
+      /* Maximum number of iteration. */
+      return LBFGSERR_MAXIMUMLINESEARCH;
+    }
+    if (*f <= ftest1 && fabs(dg) <= param->gtol * (-dginit)) {
+      /* The sufficient decrease condition and the directional derivative condition hold. */
+      return count;
+    }
+	
+    /*
+      In the first stage we seek a step for which the modified
+      function has a nonpositive value and nonnegative derivative.
+    */
+    if (stage1 && *f <= ftest1 && min2(param->ftol, param->gtol) * dginit <= dg) {
+      stage1 = 0;
+    }
+
+    /*
+      A modified function is used to predict the step only if
+      we have not obtained a step for which the modified
+      function has a nonpositive function value and nonnegative
+      derivative, and if a lower function value has been
+      obtained but the decrease is not sufficient.
+    */
+    if (stage1 && ftest1 < *f && *f <= fx) {
+      /* Define the modified function and derivative values. */
+      fm = *f - *stp * dgtest;
+      fxm = fx - stx * dgtest;
+      fym = fy - sty * dgtest;
+      dgm = dg - dgtest;
+      dgxm = dgx - dgtest;
+      dgym = dgy - dgtest;
+
+      /*
+	Call update_trial_interval() to update the interval of
+	uncertainty and to compute the new step.
+      */
+      uinfo = update_trial_interval(
+				    &stx, &fxm, &dgxm,
+				    &sty, &fym, &dgym,
+				    stp, &fm, &dgm,
+				    stmin, stmax, &brackt
+				    );
+
+      /* Reset the function and gradient values for f. */
+      fx = fxm + stx * dgtest;
+      fy = fym + sty * dgtest;
+      dgx = dgxm + dgtest;
+      dgy = dgym + dgtest;
+    } else {
+      /*
+	Call update_trial_interval() to update the interval of
+	uncertainty and to compute the new step.
+      */
+      uinfo = update_trial_interval(
+				    &stx, &fx, &dgx,
+				    &sty, &fy, &dgy,
+				    stp, f, &dg,
+				    stmin, stmax, &brackt
+				    );
+    }
+
+    /*
+      Force a sufficient decrease in the interval of uncertainty.
+    */
+    if (brackt) {
+      if (0.66 * prev_width <= fabs(sty - stx)) {
+	*stp = stx + 0.5 * (sty - stx);
+      }
+      prev_width = width;
+      width = fabs(sty - stx);
+    }
+  }
+
+  return LBFGSERR_LOGICERROR;
+}
+
+
+
+/**
+ * Define the local variables for computing minimizers.
+ */
+#define USES_MINIMIZER					\
+  lbfgsfloatval_t a, d, gamma, theta, p, q, r, s;
+
+/**
+ * Find a minimizer of an interpolated cubic function.
+ *  @param  cm      The minimizer of the interpolated cubic.
+ *  @param  u       The value of one point, u.
+ *  @param  fu      The value of f(u).
+ *  @param  du      The value of f'(u).
+ *  @param  v       The value of another point, v.
+ *  @param  fv      The value of f(v).
+ *  @param  du      The value of f'(v).
+ */
+#define CUBIC_MINIMIZER(cm, u, fu, du, v, fv, dv)	\
+  d = (v) - (u);					\
+  theta = ((fu) - (fv)) * 3 / d + (du) + (dv);		\
+  p = fabs(theta);					\
+  q = fabs(du);						\
+  r = fabs(dv);						\
+  s = max3(p, q, r);					\
+  /* gamma = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */	\
+  a = theta / s;					\
+  gamma = s * sqrt(a * a - ((du) / s) * ((dv) / s));	\
+  if ((v) < (u)) gamma = -gamma;			\
+  p = gamma - (du) + theta;				\
+  q = gamma - (du) + gamma + (dv);			\
+  r = p / q;						\
+  (cm) = (u) + r * d;
+
+/**
+ * Find a minimizer of an interpolated cubic function.
+ *  @param  cm      The minimizer of the interpolated cubic.
+ *  @param  u       The value of one point, u.
+ *  @param  fu      The value of f(u).
+ *  @param  du      The value of f'(u).
+ *  @param  v       The value of another point, v.
+ *  @param  fv      The value of f(v).
+ *  @param  du      The value of f'(v).
+ *  @param  xmin    The maximum value.
+ *  @param  xmin    The minimum value.
+ */
+#define CUBIC_MINIMIZER2(cm, u, fu, du, v, fv, dv, xmin, xmax)	\
+  d = (v) - (u);						\
+  theta = ((fu) - (fv)) * 3 / d + (du) + (dv);			\
+  p = fabs(theta);						\
+  q = fabs(du);							\
+  r = fabs(dv);							\
+  s = max3(p, q, r);						\
+  /* gamma = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */		\
+  a = theta / s;						\
+  gamma = s * sqrt(max2(0., a * a - ((du) / s) * ((dv) / s)));	\
+  if ((u) < (v)) gamma = -gamma;				\
+  p = gamma - (dv) + theta;					\
+  q = gamma - (dv) + gamma + (du);				\
+  r = p / q;							\
+  if (r < 0. && gamma != 0.) {					\
+    (cm) = (v) - r * d;						\
+  } else if (a < 0) {						\
+    (cm) = (xmax);						\
+  } else {							\
+    (cm) = (xmin);						\
+  }
+
+/**
+ * Find a minimizer of an interpolated quadratic function.
+ *  @param  qm      The minimizer of the interpolated quadratic.
+ *  @param  u       The value of one point, u.
+ *  @param  fu      The value of f(u).
+ *  @param  du      The value of f'(u).
+ *  @param  v       The value of another point, v.
+ *  @param  fv      The value of f(v).
+ */
+#define QUARD_MINIMIZER(qm, u, fu, du, v, fv)			\
+  a = (v) - (u);						\
+  (qm) = (u) + (du) / (((fu) - (fv)) / a + (du)) / 2 * a;
+
+/**
+ * Find a minimizer of an interpolated quadratic function.
+ *  @param  qm      The minimizer of the interpolated quadratic.
+ *  @param  u       The value of one point, u.
+ *  @param  du      The value of f'(u).
+ *  @param  v       The value of another point, v.
+ *  @param  dv      The value of f'(v).
+ */
+#define QUARD_MINIMIZER2(qm, u, du, v, dv)	\
+  a = (u) - (v);				\
+  (qm) = (v) + (dv) / ((dv) - (du)) * a;
+
+/**
+ * Update a safeguarded trial value and interval for line search.
+ *
+ *  The parameter x represents the step with the least function value.
+ *  The parameter t represents the current step. This function assumes
+ *  that the derivative at the point of x in the direction of the step.
+ *  If the bracket is set to true, the minimizer has been bracketed in
+ *  an interval of uncertainty with endpoints between x and y.
+ *
+ *  @param  x       The pointer to the value of one endpoint.
+ *  @param  fx      The pointer to the value of f(x).
+ *  @param  dx      The pointer to the value of f'(x).
+ *  @param  y       The pointer to the value of another endpoint.
+ *  @param  fy      The pointer to the value of f(y).
+ *  @param  dy      The pointer to the value of f'(y).
+ *  @param  t       The pointer to the value of the trial value, t.
+ *  @param  ft      The pointer to the value of f(t).
+ *  @param  dt      The pointer to the value of f'(t).
+ *  @param  tmin    The minimum value for the trial value, t.
+ *  @param  tmax    The maximum value for the trial value, t.
+ *  @param  brackt  The pointer to the predicate if the trial value is
+ *                  bracketed.
+ *  @retval int     Status value. Zero indicates a normal termination.
+ *  
+ *  @see
+ *      Jorge J. More and David J. Thuente. Line search algorithm with
+ *      guaranteed sufficient decrease. ACM Transactions on Mathematical
+ *      Software (TOMS), Vol 20, No 3, pp. 286-307, 1994.
+ */
+static int update_trial_interval(
+				 lbfgsfloatval_t *x,
+				 lbfgsfloatval_t *fx,
+				 lbfgsfloatval_t *dx,
+				 lbfgsfloatval_t *y,
+				 lbfgsfloatval_t *fy,
+				 lbfgsfloatval_t *dy,
+				 lbfgsfloatval_t *t,
+				 lbfgsfloatval_t *ft,
+				 lbfgsfloatval_t *dt,
+				 const lbfgsfloatval_t tmin,
+				 const lbfgsfloatval_t tmax,
+				 int *brackt
+				 )
+{
+  int bound;
+  int dsign = fsigndiff(dt, dx);
+  lbfgsfloatval_t mc; /* minimizer of an interpolated cubic. */
+  lbfgsfloatval_t mq; /* minimizer of an interpolated quadratic. */
+  lbfgsfloatval_t newt;   /* new trial value. */
+  USES_MINIMIZER;     /* for CUBIC_MINIMIZER and QUARD_MINIMIZER. */
+
+  /* Check the input parameters for errors. */
+  if (*brackt) {
+    if (*t <= min2(*x, *y) || max2(*x, *y) <= *t) {
+      /* The trival value t is out of the interval. */
+      return LBFGSERR_OUTOFINTERVAL;
+    }
+    if (0. <= *dx * (*t - *x)) {
+      /* The function must decrease from x. */
+      return LBFGSERR_INCREASEGRADIENT;
+    }
+    if (tmax < tmin) {
+      /* Incorrect tmin and tmax specified. */
+      return LBFGSERR_INCORRECT_TMINMAX;
+    }
+  }
+
+  /*
+    Trial value selection.
+  */
+  if (*fx < *ft) {
+    /*
+      Case 1: a higher function value.
+      The minimum is brackt. If the cubic minimizer is closer
+      to x than the quadratic one, the cubic one is taken, else
+      the average of the minimizers is taken.
+    */
+    *brackt = 1;
+    bound = 1;
+    CUBIC_MINIMIZER(mc, *x, *fx, *dx, *t, *ft, *dt);
+    QUARD_MINIMIZER(mq, *x, *fx, *dx, *t, *ft);
+    if (fabs(mc - *x) < fabs(mq - *x)) {
+      newt = mc;
+    } else {
+      newt = mc + 0.5 * (mq - mc);
+    }
+  } else if (dsign) {
+    /*
+      Case 2: a lower function value and derivatives of
+      opposite sign. The minimum is brackt. If the cubic
+      minimizer is closer to x than the quadratic (secant) one,
+      the cubic one is taken, else the quadratic one is taken.
+    */
+    *brackt = 1;
+    bound = 0;
+    CUBIC_MINIMIZER(mc, *x, *fx, *dx, *t, *ft, *dt);
+    QUARD_MINIMIZER2(mq, *x, *dx, *t, *dt);
+    if (fabs(mc - *t) > fabs(mq - *t)) {
+      newt = mc;
+    } else {
+      newt = mq;
+    }
+  } else if (fabs(*dt) < fabs(*dx)) {
+    /*
+      Case 3: a lower function value, derivatives of the
+      same sign, and the magnitude of the derivative decreases.
+      The cubic minimizer is only used if the cubic tends to
+      infinity in the direction of the minimizer or if the minimum
+      of the cubic is beyond t. Otherwise the cubic minimizer is
+      defined to be either tmin or tmax. The quadratic (secant)
+      minimizer is also computed and if the minimum is brackt
+      then the the minimizer closest to x is taken, else the one
+      farthest away is taken.
+    */
+    bound = 1;
+    CUBIC_MINIMIZER2(mc, *x, *fx, *dx, *t, *ft, *dt, tmin, tmax);
+    QUARD_MINIMIZER2(mq, *x, *dx, *t, *dt);
+    if (*brackt) {
+      if (fabs(*t - mc) < fabs(*t - mq)) {
+	newt = mc;
+      } else {
+	newt = mq;
+      }
+    } else {
+      if (fabs(*t - mc) > fabs(*t - mq)) {
+	newt = mc;
+      } else {
+	newt = mq;
+      }
+    }
+  } else {
+    /*
+      Case 4: a lower function value, derivatives of the
+      same sign, and the magnitude of the derivative does
+      not decrease. If the minimum is not brackt, the step
+      is either tmin or tmax, else the cubic minimizer is taken.
+    */
+    bound = 0;
+    if (*brackt) {
+      CUBIC_MINIMIZER(newt, *t, *ft, *dt, *y, *fy, *dy);
+    } else if (*x < *t) {
+      newt = tmax;
+    } else {
+      newt = tmin;
+    }
+  }
+
+  /*
+    Update the interval of uncertainty. This update does not
+    depend on the new step or the case analysis above.
+
+    - Case a: if f(x) < f(t),
+    x <- x, y <- t.
+    - Case b: if f(t) <= f(x) && f'(t)*f'(x) > 0,
+    x <- t, y <- y.
+    - Case c: if f(t) <= f(x) && f'(t)*f'(x) < 0, 
+    x <- t, y <- x.
+  */
+  if (*fx < *ft) {
+    /* Case a */
+    *y = *t;
+    *fy = *ft;
+    *dy = *dt;
+  } else {
+    /* Case c */
+    if (dsign) {
+      *y = *x;
+      *fy = *fx;
+      *dy = *dx;
+    }
+    /* Cases b and c */
+    *x = *t;
+    *fx = *ft;
+    *dx = *dt;
+  }
+
+  /* Clip the new trial value in [tmin, tmax]. */
+  if (tmax < newt) newt = tmax;
+  if (newt < tmin) newt = tmin;
+
+  /*
+    Redefine the new trial value if it is close to the upper bound
+    of the interval.
+  */
+  if (*brackt && bound) {
+    mq = *x + 0.66 * (*y - *x);
+    if (*x < *y) {
+      if (mq < newt) newt = mq;
+    } else {
+      if (newt < mq) newt = mq;
+    }
+  }
+
+  /* Return the new trial value. */
+  *t = newt;
+  return 0;
+}
+
+
+
+
+
+static lbfgsfloatval_t owlqn_x1norm(
+				    const lbfgsfloatval_t* x,
+				    const int start,
+				    const int n
+				    )
+{
+  int i;
+  lbfgsfloatval_t norm = 0.;
+
+  for (i = start;i < n;++i) {
+    norm += fabs(x[i]);
+  }
+
+  return norm;
+}
+
+static void owlqn_pseudo_gradient(
+				  lbfgsfloatval_t* pg,
+				  const lbfgsfloatval_t* x,
+				  const lbfgsfloatval_t* g,
+				  const int n,
+				  const lbfgsfloatval_t c,
+				  const int start,
+				  const int end
+				  )
+{
+  int i;
+
+  /* Compute the negative of gradients. */
+  for (i = 0;i < start;++i) {
+    pg[i] = g[i];
+  }
+
+  /* Compute the psuedo-gradients. */
+  for (i = start;i < end;++i) {
+    if (x[i] < 0.) {
+      /* Differentiable. */
+      pg[i] = g[i] - c;
+    } else if (0. < x[i]) {
+      /* Differentiable. */
+      pg[i] = g[i] + c;
+    } else {
+      if (g[i] < -c) {
+	/* Take the right partial derivative. */
+	pg[i] = g[i] + c;
+      } else if (c < g[i]) {
+	/* Take the left partial derivative. */
+	pg[i] = g[i] - c;
+      } else {
+	pg[i] = 0.;
+      }
+    }
+  }
+
+  for (i = end;i < n;++i) {
+    pg[i] = g[i];
+  }
+}
+
+static void owlqn_project(
+			  lbfgsfloatval_t* d,
+			  const lbfgsfloatval_t* sign,
+			  const int start,
+			  const int end
+			  )
+{
+  int i;
+
+  for (i = start;i < end;++i) {
+    if (d[i] * sign[i] <= 0) {
+      d[i] = 0;
+    }
+  }
+}
diff --git a/methods/lbfgs/lbfgs.h b/methods/lbfgs/lbfgs.h
new file mode 100644
index 0000000..7be84a7
--- /dev/null
+++ b/methods/lbfgs/lbfgs.h
@@ -0,0 +1,385 @@
+#ifndef __LBFGS_H__
+#define __LBFGS_H__
+
+#ifndef LBFGS_FLOAT
+#define LBFGS_FLOAT     64
+#endif/*LBFGS_FLOAT*/
+
+/*
+ * Activate optimization routines for IEEE754 floating point values.
+ */
+#ifndef LBFGS_IEEE_FLOAT
+#define LBFGS_IEEE_FLOAT    1
+#endif/*LBFGS_IEEE_FLOAT*/
+
+#if     LBFGS_FLOAT == 32
+typedef float lbfgsfloatval_t;
+
+#elif   LBFGS_FLOAT == 64
+typedef double lbfgsfloatval_t;
+
+#else
+#error "libLBFGS supports single (float; LBFGS_FLOAT = 32) or double (double; LBFGS_FLOAT=64) precision only."
+
+#endif
+
+
+/** 
+ * \addtogroup liblbfgs_api libLBFGS API
+ * @{
+ *
+ *  The libLBFGS API.
+ */
+
+/**
+ * Return values of lbfgs().
+ * 
+ *  Roughly speaking, a negative value indicates an error.
+ */
+enum {
+  /** L-BFGS reaches convergence. */
+  LBFGS_SUCCESS = 0,
+  LBFGS_CONVERGENCE = 0,
+  LBFGS_STOP,
+  /** The initial variables already minimize the objective function. */
+  LBFGS_ALREADY_MINIMIZED,
+  
+  /** Unknown error. */
+  LBFGSERR_UNKNOWNERROR = -1024,
+  /** Logic error. */
+  LBFGSERR_LOGICERROR,
+  /** Insufficient memory. */
+  LBFGSERR_OUTOFMEMORY,
+  /** The minimization process has been canceled. */
+  LBFGSERR_CANCELED,
+  /** Invalid number of variables specified. */
+  LBFGSERR_INVALID_N,
+  /** Invalid number of variables (for SSE) specified. */
+  LBFGSERR_INVALID_N_SSE,
+  /** The array x must be aligned to 16 (for SSE). */
+  LBFGSERR_INVALID_X_SSE,
+  /** Invalid parameter lbfgs_parameter_t::epsilon specified. */
+  LBFGSERR_INVALID_EPSILON,
+  /** Invalid parameter lbfgs_parameter_t::past specified. */
+  LBFGSERR_INVALID_TESTPERIOD,
+  /** Invalid parameter lbfgs_parameter_t::delta specified. */
+  LBFGSERR_INVALID_DELTA,
+  /** Invalid parameter lbfgs_parameter_t::linesearch specified. */
+  LBFGSERR_INVALID_LINESEARCH,
+  /** Invalid parameter lbfgs_parameter_t::max_step specified. */
+  LBFGSERR_INVALID_MINSTEP,
+  /** Invalid parameter lbfgs_parameter_t::max_step specified. */
+  LBFGSERR_INVALID_MAXSTEP,
+  /** Invalid parameter lbfgs_parameter_t::ftol specified. */
+  LBFGSERR_INVALID_FTOL,
+  /** Invalid parameter lbfgs_parameter_t::wolfe specified. */
+  LBFGSERR_INVALID_WOLFE,
+  /** Invalid parameter lbfgs_parameter_t::gtol specified. */
+  LBFGSERR_INVALID_GTOL,
+  /** Invalid parameter lbfgs_parameter_t::xtol specified. */
+  LBFGSERR_INVALID_XTOL,
+  /** Invalid parameter lbfgs_parameter_t::max_linesearch specified. */
+  LBFGSERR_INVALID_MAXLINESEARCH,
+  /** Invalid parameter lbfgs_parameter_t::orthantwise_c specified. */
+  LBFGSERR_INVALID_ORTHANTWISE,
+  /** Invalid parameter lbfgs_parameter_t::orthantwise_start specified. */
+  LBFGSERR_INVALID_ORTHANTWISE_START,
+  /** Invalid parameter lbfgs_parameter_t::orthantwise_end specified. */
+  LBFGSERR_INVALID_ORTHANTWISE_END,
+  /** The line-search step went out of the interval of uncertainty. */
+  LBFGSERR_OUTOFINTERVAL,
+  /** A logic error occurred; alternatively, the interval of uncertainty
+      became too small. */
+  LBFGSERR_INCORRECT_TMINMAX,
+  /** A rounding error occurred; alternatively, no line-search step
+      satisfies the sufficient decrease and curvature conditions. */
+  LBFGSERR_ROUNDING_ERROR,
+  /** The line-search step became smaller than lbfgs_parameter_t::min_step. */
+  LBFGSERR_MINIMUMSTEP,
+  /** The line-search step became larger than lbfgs_parameter_t::max_step. */
+  LBFGSERR_MAXIMUMSTEP,
+  /** The line-search routine reaches the maximum number of evaluations. */
+  LBFGSERR_MAXIMUMLINESEARCH,
+  /** The algorithm routine reaches the maximum number of iterations. */
+  LBFGSERR_MAXIMUMITERATION,
+  /** Relative width of the interval of uncertainty is at most
+      lbfgs_parameter_t::xtol. */
+  LBFGSERR_WIDTHTOOSMALL,
+  /** A logic error (negative line-search step) occurred. */
+  LBFGSERR_INVALIDPARAMETERS,
+  /** The current search direction increases the objective function value. */
+  LBFGSERR_INCREASEGRADIENT,
+};
+
+/**
+ * Line search algorithms.
+ */
+enum {
+  /** The default algorithm (MoreThuente method). */
+  LBFGS_LINESEARCH_DEFAULT = 0,
+  /** MoreThuente method proposd by More and Thuente. */
+  LBFGS_LINESEARCH_MORETHUENTE = 0,
+  /**
+   * Backtracking method with the Armijo condition.
+   *  The backtracking method finds the step length such that it satisfies
+   *  the sufficient decrease (Armijo) condition,
+   *    - f(x + a * d) <= f(x) + lbfgs_parameter_t::ftol * a * g(x)^T d,
+   *
+   *  where x is the current point, d is the current search direction, and
+   *  a is the step length.
+   */
+  LBFGS_LINESEARCH_BACKTRACKING_ARMIJO = 1,
+  /** The backtracking method with the defualt (regular Wolfe) condition. */
+  LBFGS_LINESEARCH_BACKTRACKING = 2,
+  /**
+   * Backtracking method with regular Wolfe condition.
+   *  The backtracking method finds the step length such that it satisfies
+   *  both the Armijo condition (LBFGS_LINESEARCH_BACKTRACKING_ARMIJO)
+   *  and the curvature condition,
+   *    - g(x + a * d)^T d >= lbfgs_parameter_t::wolfe * g(x)^T d,
+   *
+   *  where x is the current point, d is the current search direction, and
+   *  a is the step length.
+   */
+  LBFGS_LINESEARCH_BACKTRACKING_WOLFE = 2,
+  /**
+   * Backtracking method with strong Wolfe condition.
+   *  The backtracking method finds the step length such that it satisfies
+   *  both the Armijo condition (LBFGS_LINESEARCH_BACKTRACKING_ARMIJO)
+   *  and the following condition,
+   *    - |g(x + a * d)^T d| <= lbfgs_parameter_t::wolfe * |g(x)^T d|,
+   *
+   *  where x is the current point, d is the current search direction, and
+   *  a is the step length.
+   */
+  LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE = 3,
+};
+
+/**
+ * L-BFGS optimization parameters.
+ *  Call lbfgs_parameter_init() function to initialize parameters to the
+ *  default values.
+ */
+typedef struct {
+  /**
+   * The number of corrections to approximate the inverse hessian matrix.
+   *  The L-BFGS routine stores the computation results of previous \ref m
+   *  iterations to approximate the inverse hessian matrix of the current
+   *  iteration. This parameter controls the size of the limited memories
+   *  (corrections). The default value is \c 6. Values less than \c 3 are
+   *  not recommended. Large values will result in excessive computing time.
+   */
+  int             m;
+  
+  /**
+   * Epsilon for convergence test.
+   *  This parameter determines the accuracy with which the solution is to
+   *  be found. A minimization terminates when
+   *      ||g|| < \ref epsilon * max(1, ||x||),
+   *  where ||.|| denotes the Euclidean (L2) norm. The default value is
+   *  \c 1e-5.
+   */
+  lbfgsfloatval_t epsilon;
+  
+  /**
+   * Distance for delta-based convergence test.
+   *  This parameter determines the distance, in iterations, to compute
+   *  the rate of decrease of the objective function. If the value of this
+   *  parameter is zero, the library does not perform the delta-based
+   *  convergence test. The default value is \c 0.
+   */
+  int             past;
+  
+  /**
+   * Delta for convergence test.
+   *  This parameter determines the minimum rate of decrease of the
+   *  objective function. The library stops iterations when the
+   *  following condition is met:
+   *      (f' - f) / f < \ref delta,
+   *  where f' is the objective value of \ref past iterations ago, and f is
+   *  the objective value of the current iteration.
+   *  The default value is \c 0.
+   */
+  lbfgsfloatval_t delta;
+  
+  /**
+   * The maximum number of iterations.
+   *  The lbfgs() function terminates an optimization process with
+   *  ::LBFGSERR_MAXIMUMITERATION status code when the iteration count
+   *  exceedes this parameter. Setting this parameter to zero continues an
+   *  optimization process until a convergence or error. The default value
+   *  is \c 0.
+   */
+  int             max_iterations;
+  
+  /**
+   * The line search algorithm.
+   *  This parameter specifies a line search algorithm to be used by the
+   *  L-BFGS routine.
+   */
+  int             linesearch;
+  
+  /**
+   * The maximum number of trials for the line search.
+   *  This parameter controls the number of function and gradients evaluations
+   *  per iteration for the line search routine. The default value is \c 20.
+   */
+  int             max_linesearch;
+  
+  /**
+   * The minimum step of the line search routine.
+   *  The default value is \c 1e-20. This value need not be modified unless
+   *  the exponents are too large for the machine being used, or unless the
+   *  problem is extremely badly scaled (in which case the exponents should
+   *  be increased).
+   */
+  lbfgsfloatval_t min_step;
+  
+  /**
+     * The maximum step of the line search.
+     *  The default value is \c 1e+20. This value need not be modified unless
+     *  the exponents are too large for the machine being used, or unless the
+     *  problem is extremely badly scaled (in which case the exponents should
+     *  be increased).
+     */
+    lbfgsfloatval_t max_step;
+
+    /**
+     * A parameter to control the accuracy of the line search routine.
+     *  The default value is \c 1e-4. This parameter should be greater
+     *  than zero and smaller than \c 0.5.
+     */
+    lbfgsfloatval_t ftol;
+
+    /**
+     * A coefficient for the Wolfe condition.
+     *  This parameter is valid only when the backtracking line-search
+     *  algorithm is used with the Wolfe condition,
+     *  ::LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE or
+     *  ::LBFGS_LINESEARCH_BACKTRACKING_WOLFE .
+     *  The default value is \c 0.9. This parameter should be greater
+     *  the \ref ftol parameter and smaller than \c 1.0.
+     */
+    lbfgsfloatval_t wolfe;
+
+    /**
+     * A parameter to control the accuracy of the line search routine.
+     *  The default value is \c 0.9. If the function and gradient
+     *  evaluations are inexpensive with respect to the cost of the
+     *  iteration (which is sometimes the case when solving very large
+     *  problems) it may be advantageous to set this parameter to a small
+     *  value. A typical small value is \c 0.1. This parameter shuold be
+     *  greater than the \ref ftol parameter (\c 1e-4) and smaller than
+     *  \c 1.0.
+     */
+    lbfgsfloatval_t gtol;
+
+    /**
+     * The machine precision for floating-point values.
+     *  This parameter must be a positive value set by a client program to
+     *  estimate the machine precision. The line search routine will terminate
+     *  with the status code (::LBFGSERR_ROUNDING_ERROR) if the relative width
+     *  of the interval of uncertainty is less than this parameter.
+     */
+    lbfgsfloatval_t xtol;
+
+    /**
+     * Coeefficient for the L1 norm of variables.
+     *  This parameter should be set to zero for standard minimization
+     *  problems. Setting this parameter to a positive value activates
+     *  Orthant-Wise Limited-memory Quasi-Newton (OWL-QN) method, which
+     *  minimizes the objective function F(x) combined with the L1 norm |x|
+     *  of the variables, {F(x) + C |x|}. This parameter is the coeefficient
+     *  for the |x|, i.e., C. As the L1 norm |x| is not differentiable at
+     *  zero, the library modifies function and gradient evaluations from
+     *  a client program suitably; a client program thus have only to return
+     *  the function value F(x) and gradients G(x) as usual. The default value
+     *  is zero.
+     */
+    lbfgsfloatval_t orthantwise_c;
+
+    /**
+     * Start index for computing L1 norm of the variables.
+     *  This parameter is valid only for OWL-QN method
+     *  (i.e., \ref orthantwise_c != 0). This parameter b (0 <= b < N)
+     *  specifies the index number from which the library computes the
+     *  L1 norm of the variables x,
+     *      |x| := |x_{b}| + |x_{b+1}| + ... + |x_{N}| .
+     *  In other words, variables x_1, ..., x_{b-1} are not used for
+     *  computing the L1 norm. Setting b (0 < b < N), one can protect
+     *  variables, x_1, ..., x_{b-1} (e.g., a bias term of logistic
+     *  regression) from being regularized. The default value is zero.
+     */
+    int             orthantwise_start;
+
+    /**
+     * End index for computing L1 norm of the variables.
+     *  This parameter is valid only for OWL-QN method
+     *  (i.e., \ref orthantwise_c != 0). This parameter e (0 < e <= N)
+     *  specifies the index number at which the library stops computing the
+     *  L1 norm of the variables x,
+     */
+    int             orthantwise_end;
+} lbfgs_parameter_t;
+
+
+/**
+ * Callback interface to provide objective function and gradient evaluations.
+ *
+ *  The lbfgs() function call this function to obtain the values of objective
+ *  function and its gradients when needed. A client program must implement
+ *  this function to evaluate the values of the objective function and its
+ *  gradients, given current values of variables.
+ *  
+ *  @param  instance    The user data sent for lbfgs() function by the client.
+ *  @param  x           The current values of variables.
+ *  @param  g           The gradient vector. The callback function must compute
+ *                      the gradient values for the current variables.
+ *  @param  n           The number of variables.
+ *  @param  step        The current step of the line search routine.
+ *  @retval lbfgsfloatval_t The value of the objective function for the current
+ *                          variables.
+ */
+typedef lbfgsfloatval_t (*lbfgs_evaluate_t)(
+    void *instance,
+    const lbfgsfloatval_t *x,
+    lbfgsfloatval_t *g,
+    const int n,
+    const lbfgsfloatval_t step
+    );
+
+/**
+ * Callback interface to receive the progress of the optimization process.
+ *
+ *  The lbfgs() function call this function for each iteration. Implementing
+ *  this function, a client program can store or display the current progress
+ *  of the optimization process.
+ *
+ *  @param  instance    The user data sent for lbfgs() function by the client.
+ *  @param  x           The current values of variables.
+ *  @param  g           The current gradient values of variables.
+ *  @param  fx          The current value of the objective function.
+ *  @param  xnorm       The Euclidean norm of the variables.
+ *  @param  gnorm       The Euclidean norm of the gradients.
+ *  @param  step        The line-search step used for this iteration.
+ *  @param  n           The number of variables.
+ *  @param  k           The iteration count.
+ *  @param  ls          The number of evaluations called for this iteration.
+ *  @retval int         Zero to continue the optimization process. Returning a
+ *                      non-zero value will cancel the optimization process.
+ */
+typedef int (*lbfgs_progress_t)(
+    void *instance,
+    const lbfgsfloatval_t *x,
+    const lbfgsfloatval_t *g,
+    const lbfgsfloatval_t fx,
+    const lbfgsfloatval_t xnorm,
+    const lbfgsfloatval_t gnorm,
+    const lbfgsfloatval_t step,
+    int n,
+    int k,
+    int ls
+    );
+
+#endif
+
diff --git a/methods/lbfgs/lbfgs.hpp b/methods/lbfgs/lbfgs.hpp
new file mode 100644
index 0000000..112d50e
--- /dev/null
+++ b/methods/lbfgs/lbfgs.hpp
@@ -0,0 +1,169 @@
+#ifndef BFGS_METHOD
+#define BFGS_METHOD
+#define OPT_HEADER
+#include <core/optimizer.hpp>
+//#include <blitz/array.h>
+#include <limits>
+#include <cstdlib>
+#include <core/opt_traits.hpp>
+#include "../linmin/linmin.hpp"
+#include <math/num_diff.hpp>
+#include <cassert>
+#include <cmath>
+#include <ctime>
+#include <vector>
+#include <algorithm>
+#include "lbfgs.h"
+#include "lbfgs.cpp"
+/*
+ *
+*/
+#include <iostream>
+using std::cerr;
+using std::endl;
+
+namespace opt_utilities
+{
+
+  template<typename rT,typename pT>
+  lbfgsfloatval_t lbfgs_adapter(
+				void *instance,
+				const lbfgsfloatval_t *x,
+				lbfgsfloatval_t *g,
+				const int n,
+				const lbfgsfloatval_t step
+				)
+  {
+    pT px;
+    resize(px,n);
+    for(int i=0;i<n;++i)
+      {
+	set_element(px,i,x[i]);
+      }
+    
+    lbfgsfloatval_t result=((func_obj<rT,pT>*)instance)->eval(px);
+    pT grad(gradient(*static_cast<func_obj<rT,pT>*>(instance),px));
+    for(int i=0;i<n;++i)
+      {
+	g[i]=get_element(grad,i);
+      }
+    return result;
+  }
+
+
+  
+  template <typename rT,typename pT>
+  class lbfgs_method
+    :public opt_method<rT,pT>
+  {
+  public:
+    typedef pT array1d_type;
+    typedef rT T;
+  private:
+    func_obj<rT,pT>* p_fo;
+    optimizer<rT,pT>* p_optimizer;
+    
+    //typedef blitz::Array<rT,2> array2d_type;
+    
+    
+  private:
+    array1d_type start_point;
+    array1d_type end_point;
+    
+  private:
+    rT threshold;
+  private:
+    rT func(const pT& x)
+    {
+      assert(p_fo!=0);
+      return p_fo->eval(x);
+    }
+
+   
+  public:
+    lbfgs_method()
+      :threshold(1e-4)
+    {}
+
+    virtual ~lbfgs_method()
+    {
+    };
+
+    lbfgs_method(const lbfgs_method<rT,pT>& rhs)
+      :p_fo(rhs.p_fo),p_optimizer(rhs.p_optimizer),
+       start_point(rhs.start_point),
+       end_point(rhs.end_point),
+       threshold(rhs.threshold)
+    {
+    }
+
+    lbfgs_method<rT,pT>& operator=(const lbfgs_method<rT,pT>& rhs)
+    {
+      threshold=rhs.threshold;
+      p_fo=rhs.p_fo;
+      p_optimizer=rhs.p_optimizer;
+      opt_eq(start_point,rhs.start_point);
+      opt_eq(end_point,rhs.end_point);
+    }
+    
+    opt_method<rT,pT>* do_clone()const
+    {
+      return new lbfgs_method<rT,pT>(*this);
+    }
+    
+    void do_set_start_point(const array1d_type& p)
+    {
+      start_point.resize(get_size(p));
+      opt_eq(start_point,p);
+      
+    }
+
+    array1d_type do_get_start_point()const
+    {
+      return start_point;
+    }
+
+    void do_set_precision(rT t)
+    {
+      threshold=t;
+    }
+
+    rT do_get_precision()const
+    {
+      return threshold;
+    }
+
+    void do_set_optimizer(optimizer<rT,pT>& o)
+    {
+      p_optimizer=&o;
+      p_fo=p_optimizer->ptr_func_obj();
+    }
+    
+    
+    
+    pT do_optimize()
+    {
+      lbfgs_parameter_t param;
+      lbfgs_parameter_init(&param);
+      param.ftol=threshold;
+      std::vector<lbfgsfloatval_t> buffer(get_size(start_point));
+      for(int i=0;i<buffer.size();++i)
+	{
+	  buffer[i]=get_element(start_point,i);
+	}
+      lbfgsfloatval_t fx;
+      lbfgs(get_size(start_point),&buffer[0],&fx,
+	    lbfgs_adapter<rT,pT>,0,p_fo,&param);
+      for(int i=0;i<buffer.size();++i)
+	{
+	  set_element(start_point,i,buffer[i]);
+	}
+      return start_point;
+    } 
+  };
+  
+}
+
+
+#endif
+//EOF