tmp update

src/CMakeLists.txt

@@ -76,3 +76,4 @@ endmacro()
 add_sample(lbfgs_sample sample.c)
 add_sample(lbfgs_sample2 sample2.cpp)
 add_sample(lbfgs_sample3 sample3.cpp)
+add_sample(lbfgs_sample4 sample4.cpp)

src/lib/lbfgs.c

@@ -268,6 +268,7 @@ int lbfgs(
     lbfgsfloatval_t *xp = NULL;
     lbfgsfloatval_t *g = NULL, *gp = NULL, *pg = NULL; // gp (p for previous) pg (p for pesudo)
     lbfgsfloatval_t *d = NULL, *w = NULL, *pf = NULL;
+    lbfgsfloatval_t *dp = NULL; // p for preconditioned
     iteration_data_t *lm = NULL, *it = NULL;
     lbfgsfloatval_t ys, yy;
     lbfgsfloatval_t xnorm, gnorm, beta;
@@ -383,6 +384,14 @@ int lbfgs(
         goto lbfgs_exit;
     }
 
+    if (cd.proc_precondition) {
+        dp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));
+        if (dp == NULL) {
+            ret = LBFGSERR_OUTOFMEMORY;
+            goto lbfgs_exit;
+        }
+    }
+
     // 初始化计算L1模的数组
     if (param.orthantwise_c != 0.) {
         /* Allocate working space for OW-LQN. */
@@ -620,7 +629,13 @@ int lbfgs(
 
         // 我们在这里提供一个预优函数的接口 返回 d = H0^-1 * d
         if (cd.proc_precondition) {
-            cd.proc_precondition(cd.instance, x, d, n);
+            if (param.orthantwise_c == 0.) {
+                cd.proc_precondition(cd.instance, x, g, d, dp, n);
+            } else {
+                cd.proc_precondition(cd.instance, x, pg, d, dp, n);
+            }
+            
+            veccpy(d, dp, n);
         }
 
         for (i = 0;i < bound;++i) {

src/lib/lbfgs.h

@@ -464,11 +464,26 @@ typedef int (*lbfgs_progress_t)(
     );
 
 
-
+/**
+ * Callback interface to implement the preconditioning process.
+ *
+ *  The lbfgs() function call this function for each iteration. Implementing
+ *  this function, a client program can preform the preconditioning 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  d           The current values of search directions.
+ *  @param  d_pre       The values of search directions being preconditioned. 
+ *                      The callback function must compute these values.
+ *  @param  n           The number of variables.
+ */
 typedef void (*lbfgs_precondition_t)(
     void *instance,
     const lbfgsfloatval_t *x,
-    lbfgsfloatval_t *d,
+    const lbfgsfloatval_t *g,
+    const lbfgsfloatval_t *d,
+    lbfgsfloatval_t *d_pre,
     const int n
     );
 

src/sample/sample4.cpp

@@ -0,0 +1,200 @@
+#include "iostream"
+#include "cmath"
+#include "ctime"
+#include "random"
+
+#include "../lib/lbfgs.h"
+
+using std::clog;
+using std::endl;
+
+#define N 100
+
+inline lbfgsfloatval_t scalernd(lbfgsfloatval_t l, lbfgsfloatval_t h)
+{
+    return (h-l)*rand()*1.0/RAND_MAX + l;
+}
+
+inline void vecrnd(lbfgsfloatval_t *a, lbfgsfloatval_t l, lbfgsfloatval_t h, int n)
+{
+    srand(time(nullptr));
+	for (size_t i = 0; i < n; i++)
+	{
+		a[i] = (h-l)*rand()*1.0/RAND_MAX + l;
+	}
+    return;
+}
+
+inline lbfgsfloatval_t max_diff(const lbfgsfloatval_t *a, const lbfgsfloatval_t *b, int n)
+{
+	lbfgsfloatval_t tmp, max = -1.0;
+	for (int i = 0; i < n; i++)
+	{
+        tmp = sqrt((a[i] - b[i])*(a[i] - b[i]));
+		max = tmp>max?tmp:max;
+	}
+	return max;
+}
+
+class TEST_FUNC
+{
+public:
+	TEST_FUNC();
+	virtual ~TEST_FUNC();
+
+	static lbfgsfloatval_t _Func(void *instance, const lbfgsfloatval_t *x, lbfgsfloatval_t *g,
+		const int n, const lbfgsfloatval_t step)
+	{
+		return reinterpret_cast<TEST_FUNC*>(instance)->Func(x, g, n, step);
+	}
+
+	lbfgsfloatval_t Func(const lbfgsfloatval_t *x, lbfgsfloatval_t *g,
+		const int n, const lbfgsfloatval_t step);
+
+	static int _Progress(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, const lbfgs_parameter_t param,
+		int n, int k, int ls)
+	{
+		return reinterpret_cast<TEST_FUNC*>(instance)->Progress(x, g, fx, xnorm, gnorm, step, param, n, k, ls);
+	}
+
+	int Progress(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t fx,
+		const lbfgsfloatval_t xnorm, const lbfgsfloatval_t gnorm, const lbfgsfloatval_t step, const lbfgs_parameter_t param,
+		int n, int k, int ls);
+
+    static void _Precondition(void *instance, const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t *d,
+        lbfgsfloatval_t *d_pre, const int n)
+	{
+		return reinterpret_cast<TEST_FUNC*>(instance)->Precondition(x, g, d, d_pre, n);
+	}
+
+	void Precondition(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t *d,
+        lbfgsfloatval_t *d_pre, const int n);
+
+	int Routine();
+private:
+    lbfgsfloatval_t **kernel;
+	lbfgsfloatval_t *x, *m, *mid, *d;
+};
+
+TEST_FUNC::TEST_FUNC()
+{
+	kernel = new lbfgsfloatval_t* [N];
+    for (size_t i = 0; i < N; i++)
+    {
+        kernel[i] =  new lbfgsfloatval_t [N];
+    }
+    
+    x = new lbfgsfloatval_t [N];
+    m = new lbfgsfloatval_t [N];
+    mid = new lbfgsfloatval_t [N];
+    d = new lbfgsfloatval_t [N];
+
+    srand(time(nullptr));
+    for (size_t i = 0; i < N-1; i++)
+    {
+        for (size_t j = i+1; j < N; j++)
+        {
+            kernel[i][j] = scalernd(0.0, 2.0);
+            kernel[j][i] = kernel[i][j];
+        }
+
+        kernel[i][i] = scalernd(10.0, 20.0);
+    }
+
+    vecrnd(m, 1.0, 2.0, N);
+
+    for (size_t i = 0; i < N; i++)
+    {
+        d[i] = 0.0;
+        for (size_t j = 0; j < N; j++)
+        {
+           d[i] += kernel[i][j] * m[j];
+        }
+    }
+}
+
+TEST_FUNC::~TEST_FUNC()
+{
+    for (size_t i = 0; i < N; i++)
+    {
+        delete[] kernel[i];
+    }
+    delete[] kernel;
+
+    delete[] x;
+    delete[] m;
+    delete[] mid;
+    delete[] d;
+}
+
+lbfgsfloatval_t TEST_FUNC::Func(const lbfgsfloatval_t *x, lbfgsfloatval_t *g,
+		const int n, const lbfgsfloatval_t step)
+{
+    for (size_t i = 0; i < N; i++)
+    {
+        mid[i] = 0.0;
+        for (size_t j = 0; j < N; j++)
+        {
+            mid[i] += kernel[i][j] * x[j];
+        }
+        mid[i] -= d[i];
+    }
+
+    lbfgsfloatval_t fx = 0;
+    for (size_t i = 0; i < N; i++)
+    {
+        fx += mid[i]*mid[i];
+    }
+    
+    for (size_t j = 0; j < N; j++)
+    {
+        g[j] = 0.0;
+        for (size_t i = 0; i < N; i++)
+        {
+            g[j] += kernel[i][j]*mid[i];
+        }
+    }
+	return fx;
+}
+
+int TEST_FUNC::Progress(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t fx,
+		const lbfgsfloatval_t xnorm, const lbfgsfloatval_t gnorm, const lbfgsfloatval_t step, const lbfgs_parameter_t param,
+		int n, int k, int ls)
+{
+	clog << "\riteration times: " << k << " fx = " << fx << " gnorm/xnorm = " << gnorm/(xnorm>1.0?xnorm:1.0);
+	return 0;
+}
+
+void TEST_FUNC::Precondition(const lbfgsfloatval_t *x, const lbfgsfloatval_t *g, const lbfgsfloatval_t *d,
+    lbfgsfloatval_t *d_pre, const int n)
+{
+    for (size_t j = 0; j < N; j++)
+    {
+        d_pre[j] = d[j]/kernel[j][j];
+    }
+    return;
+}
+
+int TEST_FUNC::Routine()
+{
+	lbfgsfloatval_t fx;
+
+	lbfgs_parameter_t self_para;
+	lbfgs_parameter_init(&self_para);
+	self_para.epsilon = 1e-7;
+    self_para.max_linesearch = 1000;
+
+	int ret = lbfgs(N, x, &fx, _Func, _Progress, this, &self_para, _Precondition);
+
+	clog << endl << "L-BFGS optimization terminated with status: " << lbfgs_strerror(ret) << endl;
+    clog << "Maximal difference = " << max_diff(x, m, N) << endl;
+	return ret;
+}
+
+int main(int argc, char const *argv[])
+{
+	TEST_FUNC test;
+	test.Routine();
+	return 0;
+}