diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 09c96de..7374d43 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -75,4 +75,5 @@ endmacro() add_sample(lbfgs_sample sample.c) add_sample(lbfgs_sample2 sample2.cpp) -add_sample(lbfgs_sample3 sample3.cpp) \ No newline at end of file +add_sample(lbfgs_sample3 sample3.cpp) +add_sample(lbfgs_sample4 sample4.cpp) \ No newline at end of file diff --git a/src/lib/lbfgs.c b/src/lib/lbfgs.c index 2bd885d..e037e8b 100644 --- a/src/lib/lbfgs.c +++ b/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) { diff --git a/src/lib/lbfgs.h b/src/lib/lbfgs.h index 18070c8..4df1ee0 100644 --- a/src/lib/lbfgs.h +++ b/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 ); diff --git a/src/sample/sample4.cpp b/src/sample/sample4.cpp new file mode 100644 index 0000000..32dcb6b --- /dev/null +++ b/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(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(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(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; +} \ No newline at end of file