From: Dima Kogan Date: Sat, 17 Nov 2012 12:25:43 -0800 Subject: using lbfgs backend from the liblfgs package, so do not need it here Forwarded: not-needed --- MANIFEST | 5 - Makefile.PL | 2 +- arithmetic_ansi.h | 133 ------ arithmetic_sse_double.h | 279 ------------ arithmetic_sse_float.h | 289 ------------- lbfgs.c | 1105 ----------------------------------------------- lbfgs.h | 508 ---------------------- t/01-parameter.t | 2 +- 8 files changed, 2 insertions(+), 2321 deletions(-) delete mode 100644 arithmetic_ansi.h delete mode 100644 arithmetic_sse_double.h delete mode 100644 arithmetic_sse_float.h delete mode 100644 lbfgs.c delete mode 100644 lbfgs.h diff --git a/MANIFEST b/MANIFEST index 046bbed..bf91503 100644 --- a/MANIFEST +++ b/MANIFEST @@ -1,7 +1,4 @@ Algorithm-LBFGS.xs -arithmetic_ansi.h -arithmetic_sse_double.h -arithmetic_sse_float.h Changes inc/Module/AutoInstall.pm inc/Module/Install.pm @@ -15,8 +12,6 @@ inc/Test/Builder.pm inc/Test/Builder/Module.pm inc/Test/More.pm inc/Test/Number/Delta.pm -lbfgs.c -lbfgs.h lib/Algorithm/LBFGS.pm LICENSE Makefile.PL diff --git a/Makefile.PL b/Makefile.PL index 1763c0b..c53bd7e 100644 --- a/Makefile.PL +++ b/Makefile.PL @@ -15,7 +15,7 @@ include 'Test::Number::Delta'; auto_install; WriteMakefile( - LIBS => ['-lm'], + LIBS => ['-llbfgs'], INC => '-I.', OBJECT => '$(O_FILES)' ); diff --git a/arithmetic_ansi.h b/arithmetic_ansi.h deleted file mode 100644 index 1458ce9..0000000 --- a/arithmetic_ansi.h +++ /dev/null @@ -1,133 +0,0 @@ -/* - * ANSI C implementation of vector operations. - * - * Copyright (c) 2007, 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 14 2007-10-25 02:04:09Z naoaki $ */ - -#include -#include - -#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 vecnorm(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n) -{ - vecdot(s, x, x, n); - *s = (lbfgsfloatval_t)sqrt(*s); -} - -inline static void vecrnorm(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n) -{ - vecnorm(s, x, n); - *s = (lbfgsfloatval_t)(1.0 / *s); -} diff --git a/arithmetic_sse_double.h b/arithmetic_sse_double.h deleted file mode 100644 index e69f7c1..0000000 --- a/arithmetic_sse_double.h +++ /dev/null @@ -1,279 +0,0 @@ -/* - * SSE2 implementation of vector oprations (64bit double). - * - * Copyright (c) 2007, 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_sse_double.h 2 2007-10-20 01:38:42Z naoaki $ */ - -#include -#include -#include - -#if 1400 <= _MSC_VER -#include -#endif - -#ifndef _MSC_VER -#include -#endif - -inline static void* vecalloc(size_t size) -{ - void *memblock = _aligned_malloc(size, 16); - if (memblock != NULL) { - memset(memblock, 0, size); - } - return memblock; -} - -inline static void vecfree(void *memblock) -{ - _aligned_free(memblock); -} - -#define fsigndiff(x, y) \ - ((_mm_movemask_pd(_mm_set_pd(*(x), *(y))) + 1) & 0x002) - -#define vecset(x, c, n) \ -{ \ - int i; \ - __m128d XMM0 = _mm_set1_pd(c); \ - for (i = 0;i < (n);i += 8) { \ - _mm_store_pd((x)+i , XMM0); \ - _mm_store_pd((x)+i+2, XMM0); \ - _mm_store_pd((x)+i+4, XMM0); \ - _mm_store_pd((x)+i+6, XMM0); \ - } \ -} - -#define veccpy(y, x, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 8) { \ - __m128d XMM0 = _mm_load_pd((x)+i ); \ - __m128d XMM1 = _mm_load_pd((x)+i+2); \ - __m128d XMM2 = _mm_load_pd((x)+i+4); \ - __m128d XMM3 = _mm_load_pd((x)+i+6); \ - _mm_store_pd((y)+i , XMM0); \ - _mm_store_pd((y)+i+2, XMM1); \ - _mm_store_pd((y)+i+4, XMM2); \ - _mm_store_pd((y)+i+6, XMM3); \ - } \ -} - -#define vecncpy(y, x, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 8) { \ - __m128d XMM0 = _mm_setzero_pd(); \ - __m128d XMM1 = _mm_setzero_pd(); \ - __m128d XMM2 = _mm_setzero_pd(); \ - __m128d XMM3 = _mm_setzero_pd(); \ - __m128d XMM4 = _mm_load_pd((x)+i ); \ - __m128d XMM5 = _mm_load_pd((x)+i+2); \ - __m128d XMM6 = _mm_load_pd((x)+i+4); \ - __m128d XMM7 = _mm_load_pd((x)+i+6); \ - XMM0 = _mm_sub_pd(XMM0, XMM4); \ - XMM1 = _mm_sub_pd(XMM1, XMM5); \ - XMM2 = _mm_sub_pd(XMM2, XMM6); \ - XMM3 = _mm_sub_pd(XMM3, XMM7); \ - _mm_store_pd((y)+i , XMM0); \ - _mm_store_pd((y)+i+2, XMM1); \ - _mm_store_pd((y)+i+4, XMM2); \ - _mm_store_pd((y)+i+6, XMM3); \ - } \ -} - -#define vecadd(y, x, c, n) \ -{ \ - int i; \ - __m128d XMM7 = _mm_set1_pd(c); \ - for (i = 0;i < (n);i += 4) { \ - __m128d XMM0 = _mm_load_pd((x)+i ); \ - __m128d XMM1 = _mm_load_pd((x)+i+2); \ - __m128d XMM2 = _mm_load_pd((y)+i ); \ - __m128d XMM3 = _mm_load_pd((y)+i+2); \ - XMM0 = _mm_mul_pd(XMM0, XMM7); \ - XMM1 = _mm_mul_pd(XMM1, XMM7); \ - XMM2 = _mm_add_pd(XMM2, XMM0); \ - XMM3 = _mm_add_pd(XMM3, XMM1); \ - _mm_store_pd((y)+i , XMM2); \ - _mm_store_pd((y)+i+2, XMM3); \ - } \ -} - -#define vecdiff(z, x, y, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 8) { \ - __m128d XMM0 = _mm_load_pd((x)+i ); \ - __m128d XMM1 = _mm_load_pd((x)+i+2); \ - __m128d XMM2 = _mm_load_pd((x)+i+4); \ - __m128d XMM3 = _mm_load_pd((x)+i+6); \ - __m128d XMM4 = _mm_load_pd((y)+i ); \ - __m128d XMM5 = _mm_load_pd((y)+i+2); \ - __m128d XMM6 = _mm_load_pd((y)+i+4); \ - __m128d XMM7 = _mm_load_pd((y)+i+6); \ - XMM0 = _mm_sub_pd(XMM0, XMM4); \ - XMM1 = _mm_sub_pd(XMM1, XMM5); \ - XMM2 = _mm_sub_pd(XMM2, XMM6); \ - XMM3 = _mm_sub_pd(XMM3, XMM7); \ - _mm_store_pd((z)+i , XMM0); \ - _mm_store_pd((z)+i+2, XMM1); \ - _mm_store_pd((z)+i+4, XMM2); \ - _mm_store_pd((z)+i+6, XMM3); \ - } \ -} - -#define vecscale(y, c, n) \ -{ \ - int i; \ - __m128d XMM7 = _mm_set1_pd(c); \ - for (i = 0;i < (n);i += 4) { \ - __m128d XMM0 = _mm_load_pd((y)+i ); \ - __m128d XMM1 = _mm_load_pd((y)+i+2); \ - XMM0 = _mm_mul_pd(XMM0, XMM7); \ - XMM1 = _mm_mul_pd(XMM1, XMM7); \ - _mm_store_pd((y)+i , XMM0); \ - _mm_store_pd((y)+i+2, XMM1); \ - } \ -} - -#define vecmul(y, x, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 8) { \ - __m128d XMM0 = _mm_load_pd((x)+i ); \ - __m128d XMM1 = _mm_load_pd((x)+i+2); \ - __m128d XMM2 = _mm_load_pd((x)+i+4); \ - __m128d XMM3 = _mm_load_pd((x)+i+6); \ - __m128d XMM4 = _mm_load_pd((y)+i ); \ - __m128d XMM5 = _mm_load_pd((y)+i+2); \ - __m128d XMM6 = _mm_load_pd((y)+i+4); \ - __m128d XMM7 = _mm_load_pd((y)+i+6); \ - XMM4 = _mm_mul_pd(XMM4, XMM0); \ - XMM5 = _mm_mul_pd(XMM5, XMM1); \ - XMM6 = _mm_mul_pd(XMM6, XMM2); \ - XMM7 = _mm_mul_pd(XMM7, XMM3); \ - _mm_store_pd((y)+i , XMM4); \ - _mm_store_pd((y)+i+2, XMM5); \ - _mm_store_pd((y)+i+4, XMM6); \ - _mm_store_pd((y)+i+6, XMM7); \ - } \ -} - - - -#if 3 <= __SSE__ -/* - Horizontal add with haddps SSE3 instruction. The work register (rw) - is unused. - */ -#define __horizontal_sum(r, rw) \ - r = _mm_hadd_ps(r, r); \ - r = _mm_hadd_ps(r, r); - -#else -/* - Horizontal add with SSE instruction. The work register (rw) is used. - */ -#define __horizontal_sum(r, rw) \ - rw = r; \ - r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(1, 0, 3, 2)); \ - r = _mm_add_ps(r, rw); \ - rw = r; \ - r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(2, 3, 0, 1)); \ - r = _mm_add_ps(r, rw); - -#endif - -#define vecdot(s, x, y, n) \ -{ \ - int i; \ - __m128d XMM0 = _mm_setzero_pd(); \ - __m128d XMM1 = _mm_setzero_pd(); \ - __m128d XMM2, XMM3, XMM4, XMM5; \ - for (i = 0;i < (n);i += 4) { \ - XMM2 = _mm_load_pd((x)+i ); \ - XMM3 = _mm_load_pd((x)+i+2); \ - XMM4 = _mm_load_pd((y)+i ); \ - XMM5 = _mm_load_pd((y)+i+2); \ - XMM2 = _mm_mul_pd(XMM2, XMM4); \ - XMM3 = _mm_mul_pd(XMM3, XMM5); \ - XMM0 = _mm_add_pd(XMM0, XMM2); \ - XMM1 = _mm_add_pd(XMM1, XMM3); \ - } \ - XMM0 = _mm_add_pd(XMM0, XMM1); \ - XMM1 = _mm_shuffle_pd(XMM0, XMM0, _MM_SHUFFLE2(1, 1)); \ - XMM0 = _mm_add_pd(XMM0, XMM1); \ - _mm_store_sd((s), XMM0); \ -} - -#define vecnorm(s, x, n) \ -{ \ - int i; \ - __m128d XMM0 = _mm_setzero_pd(); \ - __m128d XMM1 = _mm_setzero_pd(); \ - __m128d XMM2, XMM3, XMM4, XMM5; \ - for (i = 0;i < (n);i += 4) { \ - XMM2 = _mm_load_pd((x)+i ); \ - XMM3 = _mm_load_pd((x)+i+2); \ - XMM4 = XMM2; \ - XMM5 = XMM3; \ - XMM2 = _mm_mul_pd(XMM2, XMM4); \ - XMM3 = _mm_mul_pd(XMM3, XMM5); \ - XMM0 = _mm_add_pd(XMM0, XMM2); \ - XMM1 = _mm_add_pd(XMM1, XMM3); \ - } \ - XMM0 = _mm_add_pd(XMM0, XMM1); \ - XMM1 = _mm_shuffle_pd(XMM0, XMM0, _MM_SHUFFLE2(1, 1)); \ - XMM0 = _mm_add_pd(XMM0, XMM1); \ - XMM0 = _mm_sqrt_pd(XMM0); \ - _mm_store_sd((s), XMM0); \ -} - - -#define vecrnorm(s, x, n) \ -{ \ - int i; \ - __m128d XMM0 = _mm_setzero_pd(); \ - __m128d XMM1 = _mm_setzero_pd(); \ - __m128d XMM2, XMM3, XMM4, XMM5; \ - for (i = 0;i < (n);i += 4) { \ - XMM2 = _mm_load_pd((x)+i ); \ - XMM3 = _mm_load_pd((x)+i+2); \ - XMM4 = XMM2; \ - XMM5 = XMM3; \ - XMM2 = _mm_mul_pd(XMM2, XMM4); \ - XMM3 = _mm_mul_pd(XMM3, XMM5); \ - XMM0 = _mm_add_pd(XMM0, XMM2); \ - XMM1 = _mm_add_pd(XMM1, XMM3); \ - } \ - XMM2 = _mm_set1_pd(1.0); \ - XMM0 = _mm_add_pd(XMM0, XMM1); \ - XMM1 = _mm_shuffle_pd(XMM0, XMM0, _MM_SHUFFLE2(1, 1)); \ - XMM0 = _mm_add_pd(XMM0, XMM1); \ - XMM0 = _mm_sqrt_pd(XMM0); \ - XMM2 = _mm_div_pd(XMM2, XMM0); \ - _mm_store_sd((s), XMM2); \ -} diff --git a/arithmetic_sse_float.h b/arithmetic_sse_float.h deleted file mode 100644 index 315d778..0000000 --- a/arithmetic_sse_float.h +++ /dev/null @@ -1,289 +0,0 @@ -/* - * SSE/SSE3 implementation of vector oprations (32bit float). - * - * Copyright (c) 2007, 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_sse_float.h 2 2007-10-20 01:38:42Z naoaki $ */ - -#include -#include -#include - -#if 1400 <= _MSC_VER -#include -#endif - -/* this block is added by laye */ -#ifndef _MSC_VER -#include -#include -#endif - -#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 = _aligned_malloc(size, 16); - if (memblock != NULL) { - memset(memblock, 0, size); - } - return memblock; -} - -inline static void vecfree(void *memblock) -{ - _aligned_free(memblock); -} - -#define vecset(x, c, n) \ -{ \ - int i; \ - __m128 XMM0 = _mm_set_ps1(c); \ - for (i = 0;i < (n);i += 16) { \ - _mm_store_ps((x)+i , XMM0); \ - _mm_store_ps((x)+i+ 4, XMM0); \ - _mm_store_ps((x)+i+ 8, XMM0); \ - _mm_store_ps((x)+i+12, XMM0); \ - } \ -} - -#define veccpy(y, x, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 16) { \ - __m128 XMM0 = _mm_load_ps((x)+i ); \ - __m128 XMM1 = _mm_load_ps((x)+i+ 4); \ - __m128 XMM2 = _mm_load_ps((x)+i+ 8); \ - __m128 XMM3 = _mm_load_ps((x)+i+12); \ - _mm_store_ps((y)+i , XMM0); \ - _mm_store_ps((y)+i+ 4, XMM1); \ - _mm_store_ps((y)+i+ 8, XMM2); \ - _mm_store_ps((y)+i+12, XMM3); \ - } \ -} - -#define vecncpy(y, x, n) \ -{ \ - int i; \ - const uint32_t mask = 0x80000000; \ - __m128 XMM4 = _mm_load_ps1((float*)&mask); \ - for (i = 0;i < (n);i += 16) { \ - __m128 XMM0 = _mm_load_ps((x)+i ); \ - __m128 XMM1 = _mm_load_ps((x)+i+ 4); \ - __m128 XMM2 = _mm_load_ps((x)+i+ 8); \ - __m128 XMM3 = _mm_load_ps((x)+i+12); \ - XMM0 = _mm_xor_ps(XMM0, XMM4); \ - XMM1 = _mm_xor_ps(XMM1, XMM4); \ - XMM2 = _mm_xor_ps(XMM2, XMM4); \ - XMM3 = _mm_xor_ps(XMM3, XMM4); \ - _mm_store_ps((y)+i , XMM0); \ - _mm_store_ps((y)+i+ 4, XMM1); \ - _mm_store_ps((y)+i+ 8, XMM2); \ - _mm_store_ps((y)+i+12, XMM3); \ - } \ -} - -#define vecadd(y, x, c, n) \ -{ \ - int i; \ - __m128 XMM7 = _mm_set_ps1(c); \ - for (i = 0;i < (n);i += 8) { \ - __m128 XMM0 = _mm_load_ps((x)+i ); \ - __m128 XMM1 = _mm_load_ps((x)+i+4); \ - __m128 XMM2 = _mm_load_ps((y)+i ); \ - __m128 XMM3 = _mm_load_ps((y)+i+4); \ - XMM0 = _mm_mul_ps(XMM0, XMM7); \ - XMM1 = _mm_mul_ps(XMM1, XMM7); \ - XMM2 = _mm_add_ps(XMM2, XMM0); \ - XMM3 = _mm_add_ps(XMM3, XMM1); \ - _mm_store_ps((y)+i , XMM2); \ - _mm_store_ps((y)+i+4, XMM3); \ - } \ -} - -#define vecdiff(z, x, y, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 16) { \ - __m128 XMM0 = _mm_load_ps((x)+i ); \ - __m128 XMM1 = _mm_load_ps((x)+i+ 4); \ - __m128 XMM2 = _mm_load_ps((x)+i+ 8); \ - __m128 XMM3 = _mm_load_ps((x)+i+12); \ - __m128 XMM4 = _mm_load_ps((y)+i ); \ - __m128 XMM5 = _mm_load_ps((y)+i+ 4); \ - __m128 XMM6 = _mm_load_ps((y)+i+ 8); \ - __m128 XMM7 = _mm_load_ps((y)+i+12); \ - XMM0 = _mm_sub_ps(XMM0, XMM4); \ - XMM1 = _mm_sub_ps(XMM1, XMM5); \ - XMM2 = _mm_sub_ps(XMM2, XMM6); \ - XMM3 = _mm_sub_ps(XMM3, XMM7); \ - _mm_store_ps((z)+i , XMM0); \ - _mm_store_ps((z)+i+ 4, XMM1); \ - _mm_store_ps((z)+i+ 8, XMM2); \ - _mm_store_ps((z)+i+12, XMM3); \ - } \ -} - -#define vecscale(y, c, n) \ -{ \ - int i; \ - __m128 XMM7 = _mm_set_ps1(c); \ - for (i = 0;i < (n);i += 8) { \ - __m128 XMM0 = _mm_load_ps((y)+i ); \ - __m128 XMM1 = _mm_load_ps((y)+i+4); \ - XMM0 = _mm_mul_ps(XMM0, XMM7); \ - XMM1 = _mm_mul_ps(XMM1, XMM7); \ - _mm_store_ps((y)+i , XMM0); \ - _mm_store_ps((y)+i+4, XMM1); \ - } \ -} - -#define vecmul(y, x, n) \ -{ \ - int i; \ - for (i = 0;i < (n);i += 16) { \ - __m128 XMM0 = _mm_load_ps((x)+i ); \ - __m128 XMM1 = _mm_load_ps((x)+i+ 4); \ - __m128 XMM2 = _mm_load_ps((x)+i+ 8); \ - __m128 XMM3 = _mm_load_ps((x)+i+12); \ - __m128 XMM4 = _mm_load_ps((y)+i ); \ - __m128 XMM5 = _mm_load_ps((y)+i+ 4); \ - __m128 XMM6 = _mm_load_ps((y)+i+ 8); \ - __m128 XMM7 = _mm_load_ps((y)+i+12); \ - XMM4 = _mm_mul_ps(XMM4, XMM0); \ - XMM5 = _mm_mul_ps(XMM5, XMM1); \ - XMM6 = _mm_mul_ps(XMM6, XMM2); \ - XMM7 = _mm_mul_ps(XMM7, XMM3); \ - _mm_store_ps((y)+i , XMM4); \ - _mm_store_ps((y)+i+ 4, XMM5); \ - _mm_store_ps((y)+i+ 8, XMM6); \ - _mm_store_ps((y)+i+12, XMM7); \ - } \ -} - - - -#if 3 <= __SSE__ -/* - Horizontal add with haddps SSE3 instruction. The work register (rw) - is unused. - */ -#define __horizontal_sum(r, rw) \ - r = _mm_hadd_ps(r, r); \ - r = _mm_hadd_ps(r, r); - -#else -/* - Horizontal add with SSE instruction. The work register (rw) is used. - */ -#define __horizontal_sum(r, rw) \ - rw = r; \ - r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(1, 0, 3, 2)); \ - r = _mm_add_ps(r, rw); \ - rw = r; \ - r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(2, 3, 0, 1)); \ - r = _mm_add_ps(r, rw); - -#endif - -#define vecdot(s, x, y, n) \ -{ \ - int i; \ - __m128 XMM0 = _mm_setzero_ps(); \ - __m128 XMM1 = _mm_setzero_ps(); \ - __m128 XMM2, XMM3, XMM4, XMM5; \ - for (i = 0;i < (n);i += 8) { \ - XMM2 = _mm_load_ps((x)+i ); \ - XMM3 = _mm_load_ps((x)+i+4); \ - XMM4 = _mm_load_ps((y)+i ); \ - XMM5 = _mm_load_ps((y)+i+4); \ - XMM2 = _mm_mul_ps(XMM2, XMM4); \ - XMM3 = _mm_mul_ps(XMM3, XMM5); \ - XMM0 = _mm_add_ps(XMM0, XMM2); \ - XMM1 = _mm_add_ps(XMM1, XMM3); \ - } \ - XMM0 = _mm_add_ps(XMM0, XMM1); \ - __horizontal_sum(XMM0, XMM1); \ - _mm_store_ss((s), XMM0); \ -} - -#define vecnorm(s, x, n) \ -{ \ - int i; \ - __m128 XMM0 = _mm_setzero_ps(); \ - __m128 XMM1 = _mm_setzero_ps(); \ - __m128 XMM2, XMM3; \ - for (i = 0;i < (n);i += 8) { \ - XMM2 = _mm_load_ps((x)+i ); \ - XMM3 = _mm_load_ps((x)+i+4); \ - XMM2 = _mm_mul_ps(XMM2, XMM2); \ - XMM3 = _mm_mul_ps(XMM3, XMM3); \ - XMM0 = _mm_add_ps(XMM0, XMM2); \ - XMM1 = _mm_add_ps(XMM1, XMM3); \ - } \ - XMM0 = _mm_add_ps(XMM0, XMM1); \ - __horizontal_sum(XMM0, XMM1); \ - XMM2 = XMM0; \ - XMM1 = _mm_rsqrt_ss(XMM0); \ - XMM3 = XMM1; \ - XMM1 = _mm_mul_ss(XMM1, XMM1); \ - XMM1 = _mm_mul_ss(XMM1, XMM3); \ - XMM1 = _mm_mul_ss(XMM1, XMM0); \ - XMM1 = _mm_mul_ss(XMM1, _mm_set_ss(-0.5f)); \ - XMM3 = _mm_mul_ss(XMM3, _mm_set_ss(1.5f)); \ - XMM3 = _mm_add_ss(XMM3, XMM1); \ - XMM3 = _mm_mul_ss(XMM3, XMM2); \ - _mm_store_ss((s), XMM3); \ -} - -#define vecrnorm(s, x, n) \ -{ \ - int i; \ - __m128 XMM0 = _mm_setzero_ps(); \ - __m128 XMM1 = _mm_setzero_ps(); \ - __m128 XMM2, XMM3; \ - for (i = 0;i < (n);i += 16) { \ - XMM2 = _mm_load_ps((x)+i ); \ - XMM3 = _mm_load_ps((x)+i+4); \ - XMM2 = _mm_mul_ps(XMM2, XMM2); \ - XMM3 = _mm_mul_ps(XMM3, XMM3); \ - XMM0 = _mm_add_ps(XMM0, XMM2); \ - XMM1 = _mm_add_ps(XMM1, XMM3); \ - } \ - XMM0 = _mm_add_ps(XMM0, XMM1); \ - __horizontal_sum(XMM0, XMM1); \ - XMM2 = XMM0; \ - XMM1 = _mm_rsqrt_ss(XMM0); \ - XMM3 = XMM1; \ - XMM1 = _mm_mul_ss(XMM1, XMM1); \ - XMM1 = _mm_mul_ss(XMM1, XMM3); \ - XMM1 = _mm_mul_ss(XMM1, XMM0); \ - XMM1 = _mm_mul_ss(XMM1, _mm_set_ss(-0.5f)); \ - XMM3 = _mm_mul_ss(XMM3, _mm_set_ss(1.5f)); \ - XMM3 = _mm_add_ss(XMM3, XMM1); \ - _mm_store_ss((s), XMM3); \ -} diff --git a/lbfgs.c b/lbfgs.c deleted file mode 100644 index f04f4e5..0000000 --- a/lbfgs.c +++ /dev/null @@ -1,1105 +0,0 @@ -/* - * Limited memory BFGS (L-BFGS). - * - * Copyright (c) 1990, Jorge Nocedal - * Copyright (c) 2007, 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: lbfgs.c 56 2007-12-16 08:01:47Z naoaki $ */ - -/* -This library is a C port of the FORTRAN implementation of Limited-memory -Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method written by Jorge Nocedal. -The original FORTRAN source code is available at: -http://www.ece.northwestern.edu/~nocedal/lbfgs.html - -The L-BFGS algorithm is described in: - - Jorge Nocedal. - Updating Quasi-Newton Matrices with Limited Storage. - Mathematics of Computation, Vol. 35, No. 151, pp. 773--782, 1980. - - Dong C. Liu and Jorge Nocedal. - On the limited memory BFGS method for large scale optimization. - Mathematical Programming B, Vol. 45, No. 3, pp. 503-528, 1989. - -The line search algorithms used in this implementation are described in: - - John E. Dennis and Robert B. Schnabel. - Numerical Methods for Unconstrained Optimization and Nonlinear - Equations, Englewood Cliffs, 1983. - - 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. - -This library also implements Orthant-Wise Limited-memory Quasi-Newton (OW-LQN) -method presented in: - - Galen Andrew and Jianfeng Gao. - Scalable training of L1-regularized log-linear models. - In Proceedings of the 24th International Conference on Machine - Learning (ICML 2007), pp. 33-40, 2007. - -I would like to thank the original author, Jorge Nocedal, who has been -distributing the effieicnt and explanatory implementation in an open source -licence. -*/ - -#ifdef HAVE_CONFIG_H -#include -#endif/*HAVE_CONFIG_H*/ - -#include -#include -#include - -#include - -#ifdef _MSC_VER -#define inline __inline -typedef unsigned int uint32_t; -#endif/*_MSC_VER*/ - -/* this block is added by laye */ -#ifndef _MSC_VER -#include -#endif/*_MSC_VER*/ - -#if defined(USE_SSE) && defined(__SSE__) && LBFGS_FLOAT == 32 -/* Use SSE optimization for 32bit float precision. */ -#include "arithmetic_sse_float.h" - -#elif defined(USE_SSE) && defined(__SSE__) && LBFGS_FLOAT == 64 -/* Use SSE2 optimization for 64bit double precision. */ -#include "arithmetic_sse_double.h" - -#else -/* No CPU specific optimization. */ -#include "arithmetic_ansi.h" - -#endif - -#define min2(a, b) ((a) <= (b) ? (a) : (b)) -#define max2(a, b) ((a) >= (b) ? (a) : (b)) -#define max3(a, b, c) max2(max2((a), (b)), (c)); - -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, 20, - 1e-20, 1e20, 1e-4, 0.9, 1.0e-16, - 0.0, -}; - -/* Forward function declarations. */ - -static int line_search_backtracking( - int n, - lbfgsfloatval_t *x, - lbfgsfloatval_t *f, - lbfgsfloatval_t *g, - lbfgsfloatval_t *s, - lbfgsfloatval_t *stp, - lbfgsfloatval_t *wa, - lbfgs_evaluate_t proc_evaluate, - void *instance, - const lbfgs_parameter_t *param - ); - -static int line_search( - int n, - lbfgsfloatval_t *x, - lbfgsfloatval_t *f, - lbfgsfloatval_t *g, - lbfgsfloatval_t *s, - lbfgsfloatval_t *stp, - lbfgsfloatval_t *wa, - lbfgs_evaluate_t proc_evaluate, - void *instance, - 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 - ); - - - - -void lbfgs_parameter_init(lbfgs_parameter_t *param) -{ - memcpy(param, &_defparam, sizeof(*param)); -} - -int lbfgs( - const 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. */ - const lbfgs_parameter_t* param = (_param != NULL) ? _param : &_defparam; - const int m = param->m; - - lbfgsfloatval_t *xp = NULL, *g = NULL, *gp = NULL, *d = NULL, *w = NULL; - iteration_data_t *lm = NULL, *it = NULL; - lbfgsfloatval_t ys, yy; - lbfgsfloatval_t norm, xnorm, gnorm, beta; - lbfgsfloatval_t fx = 0.; - - /* Check the input parameters for errors. */ - if (n <= 0) { - return LBFGSERR_INVALID_N; - } -#if defined(USE_SSE) && defined(__SSE__) - if (n % 8 != 0) { - return LBFGSERR_INVALID_N_SSE; - } -#endif/*defined(__SSE__)*/ - 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->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; - } - - /* 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; - } - - /* 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; - } - } - - /* Evaluate the function value and its gradient. */ - fx = proc_evaluate(instance, x, g, n, 0); - if (0. < param->orthantwise_c) { - /* Compute L1-regularization factor and add it to the object value. */ - norm = 0.; - for (i = 0;i < n;++i) { - norm += fabs(x[i]); - } - fx += norm * param->orthantwise_c; - } - - /* We assume the initial hessian matrix H_0 as the identity matrix. */ - if (param->orthantwise_c == 0.) { - vecncpy(d, g, n); - } else { - /* Compute the negative of psuedo-gradients. */ - for (i = 0;i < n;++i) { - if (x[i] < 0.) { - /* Differentiable. */ - d[i] = -g[i] + param->orthantwise_c; - } else if (0. < x[i]) { - /* Differentiable. */ - d[i] = -g[i] - param->orthantwise_c; - } else { - if (g[i] < -param->orthantwise_c) { - /* Take the right partial derivative. */ - d[i] = -g[i] - param->orthantwise_c; - } else if (param->orthantwise_c < g[i]) { - /* Take the left partial derivative. */ - d[i] = -g[i] + param->orthantwise_c; - } else { - d[i] = 0.; - } - } - } - } - - /* Compute the initial step: - step = 1.0 / sqrt(vecdot(d, d, n)) - */ - vecrnorm(&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. */ - ls = line_search( - n, x, &fx, g, d, &step, w, proc_evaluate, instance, param); - if (ls < 0) { - ret = ls; - goto lbfgs_exit; - } - - /* Compute x and g norms. */ - vecnorm(&gnorm, g, n); - vecnorm(&xnorm, x, n); - - /* Report the progress. */ - if (proc_progress) { - if (ret = proc_progress(instance, x, g, fx, xnorm, gnorm, step, 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 = 0; - break; - } - - 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; - - if (param->orthantwise_c == 0.) { - /* Compute the negative of gradients. */ - vecncpy(d, g, n); - } else { - /* Compute the negative of psuedo-gradients. */ - for (i = 0;i < n;++i) { - if (x[i] < 0.) { - /* Differentiable. */ - d[i] = -g[i] + param->orthantwise_c; - } else if (0. < x[i]) { - /* Differentiable. */ - d[i] = -g[i] - param->orthantwise_c; - } else { - if (g[i] < -param->orthantwise_c) { - /* Take the right partial derivative. */ - d[i] = -g[i] - param->orthantwise_c; - } else if (param->orthantwise_c < g[i]) { - /* Take the left partial derivative. */ - d[i] = -g[i] + param->orthantwise_c; - } else { - d[i] = 0.; - } - } - } - /* Store the steepest direction.*/ - veccpy(w, d, 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 = 0;i < n;++i) { - if (d[i] * w[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; - } - - /* 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(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, - lbfgsfloatval_t *xp, - lbfgs_evaluate_t proc_evaluate, - void *instance, - const lbfgs_parameter_t *param - ) -{ - int i, ret = 0, count = 0; - lbfgsfloatval_t width = 0.5, norm = 0.; - lbfgsfloatval_t finit, dginit = 0., dg, dgtest; - - /* Check the input parameters for errors. */ - if (*stp <= 0.) { - return LBFGSERR_INVALIDPARAMETERS; - } - - /* Compute the initial gradient in the search direction. */ - if (param->orthantwise_c != 0.) { - /* Use psuedo-gradients for orthant-wise updates. */ - for (i = 0;i < n;++i) { - /* Notice that: - (-s[i] < 0) <==> (g[i] < -param->orthantwise_c) - (-s[i] > 0) <==> (param->orthantwise_c < g[i]) - as the result of the lbfgs() function for orthant-wise updates. - */ - if (s[i] != 0.) { - if (x[i] < 0.) { - /* Differentiable. */ - dginit += s[i] * (g[i] - param->orthantwise_c); - } else if (0. < x[i]) { - /* Differentiable. */ - dginit += s[i] * (g[i] + param->orthantwise_c); - } else if (s[i] < 0.) { - /* Take the left partial derivative. */ - dginit += s[i] * (g[i] - param->orthantwise_c); - } else if (0. < s[i]) { - /* Take the right partial derivative. */ - dginit += s[i] * (g[i] + param->orthantwise_c); - } - } - } - } else { - 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; - - /* Copy the value of x to the work area. */ - veccpy(xp, x, n); - - for (;;) { - veccpy(x, xp, n); - vecadd(x, s, *stp, n); - - if (param->orthantwise_c != 0.) { - /* The current point is projected onto the orthant of the initial one. */ - for (i = 0;i < n;++i) { - if (x[i] * xp[i] < 0.) { - x[i] = 0.; - } - } - } - - /* Evaluate the function and gradient values. */ - *f = proc_evaluate(instance, x, g, n, *stp); - if (0. < param->orthantwise_c) { - /* Compute L1-regularization factor and add it to the object value. */ - norm = 0.; - for (i = 0;i < n;++i) { - norm += fabs(x[i]); - } - *f += norm * param->orthantwise_c; - } - - vecdot(&dg, g, s, n); - ++count; - - if (*f <= finit + *stp * dgtest) { - /* The sufficient decrease condition. */ - return count; - } - if (*stp < param->min_step) { - /* The step is the minimum value. */ - ret = LBFGSERR_MINIMUMSTEP; - break; - } - if (param->max_linesearch <= count) { - /* Maximum number of iteration. */ - ret = LBFGSERR_MAXIMUMLINESEARCH; - break; - } - - (*stp) *= width; - } - - /* Revert to the previous position. */ - veccpy(x, xp, n); - return ret; -} - - - -static int line_search( - int n, - lbfgsfloatval_t *x, - lbfgsfloatval_t *f, - lbfgsfloatval_t *g, - lbfgsfloatval_t *s, - lbfgsfloatval_t *stp, - lbfgsfloatval_t *wa, - lbfgs_evaluate_t proc_evaluate, - void *instance, - const lbfgs_parameter_t *param - ) -{ - int i, count = 0; - int brackt, stage1, uinfo = 0; - lbfgsfloatval_t dg, norm; - 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. */ - if (param->orthantwise_c != 0.) { - /* Use psuedo-gradients for orthant-wise updates. */ - dginit = 0.; - for (i = 0;i < n;++i) { - /* Notice that: - (-s[i] < 0) <==> (g[i] < -param->orthantwise_c) - (-s[i] > 0) <==> (param->orthantwise_c < g[i]) - as the result of the lbfgs() function for orthant-wise updates. - */ - if (s[i] != 0.) { - if (x[i] < 0.) { - /* Differentiable. */ - dginit += s[i] * (g[i] - param->orthantwise_c); - } else if (0. < x[i]) { - /* Differentiable. */ - dginit += s[i] * (g[i] + param->orthantwise_c); - } else if (s[i] < 0.) { - /* Take the left partial derivative. */ - dginit += s[i] * (g[i] - param->orthantwise_c); - } else if (0. < s[i]) { - /* Take the right partial derivative. */ - dginit += s[i] * (g[i] + param->orthantwise_c); - } - } - } - } else { - 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; - - /* Copy the value of x to the work area. */ - veccpy(wa, x, n); - - /* - 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, wa, n); - vecadd(x, s, *stp, n); - - if (param->orthantwise_c != 0.) { - /* The current point is projected onto the orthant of the previous one. */ - for (i = 0;i < n;++i) { - if (x[i] * wa[i] < 0.) { - x[i] = 0.; - } - } - } - - /* Evaluate the function and gradient values. */ - *f = proc_evaluate(instance, x, g, n, *stp); - if (0. < param->orthantwise_c) { - /* Compute L1-regularization factor and add it to the object value. */ - norm = 0.; - for (i = 0;i < n;++i) { - norm += fabs(x[i]); - } - *f += norm * param->orthantwise_c; - } - - ++count; - - vecdot(&dg, g, s, n); - ftest1 = finit + *stp * dgtest; - - /* 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; -} diff --git a/lbfgs.h b/lbfgs.h deleted file mode 100644 index 75ef571..0000000 --- a/lbfgs.h +++ /dev/null @@ -1,508 +0,0 @@ -/* - * C library of Limited memory BFGS (L-BFGS). - * - * Copyright (c) 1990, Jorge Nocedal - * Copyright (c) 2007, 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: lbfgs.h 58 2007-12-16 09:25:33Z naoaki $ */ - -#ifndef __LBFGS_H__ -#define __LBFGS_H__ - -#ifdef __cplusplus -extern "C" { -#endif/*__cplusplus*/ - -/* - * The default precision of floating point values is 64bit (double). - */ -#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(). - */ -enum { - /** False value. */ - LBFGSFALSE = 0, - /** True value. */ - LBFGSTRUE, - - /** 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, - /** 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::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, - /** 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, -}; - -/** - * 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; - - /** - * 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 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 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 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 - * modify 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; -} 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 - ); - -/* -A user must implement a function compatible with ::lbfgs_evaluate_t (evaluation -callback) and pass the pointer to the callback function to lbfgs() arguments. -Similarly, a user can implement a function compatible with ::lbfgs_progress_t -(progress callback) to obtain the current progress (e.g., variables, function -value, ||G||, etc) and to cancel the iteration process if necessary. -Implementation of a progress callback is optional: a user can pass \c NULL if -progress notification is not necessary. - -In addition, a user must preserve two requirements: - - The number of variables must be multiples of 16 (this is not 4). - - The memory block of variable array ::x must be aligned to 16. - -This algorithm terminates an optimization -when: - - ||G|| < \epsilon \cdot \max(1, ||x||) . - -In this formula, ||.|| denotes the Euclidean norm. -*/ - -/** - * Start a L-BFGS optimization. - * - * @param n The number of variables. - * @param x The array of variables. A client program can set - * default values for the optimization and receive the - * optimization result through this array. - * @param ptr_fx The pointer to the variable that receives the final - * value of the objective function for the variables. - * This argument can be set to \c NULL if the final - * value of the objective function is unnecessary. - * @param proc_evaluate The callback function to provide function and - * gradient evaluations given a current values of - * variables. A client program must implement a - * callback function compatible with \ref - * lbfgs_evaluate_t and pass the pointer to the - * callback function. - * @param proc_progress The callback function to receive the progress - * (the number of iterations, the current value of - * the objective function) of the minimization - * process. This argument can be set to \c NULL if - * a progress report is unnecessary. - * @param instance A user data for the client program. The callback - * functions will receive the value of this argument. - * @param param The pointer to a structure representing parameters for - * L-BFGS optimization. A client program can set this - * parameter to \c NULL to use the default parameters. - * Call lbfgs_parameter_init() function to fill a - * structure with the default values. - * @retval int The status code. This function returns zero if the - * minimization process terminates without an error. A - * non-zero value indicates an error. - */ -int lbfgs( - const 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 - ); - -/** - * Initialize L-BFGS parameters to the default values. - * - * Call this function to fill a parameter structure with the default values - * and overwrite parameter values if necessary. - * - * @param param The pointer to the parameter structure. - */ -void lbfgs_parameter_init(lbfgs_parameter_t *param); - -/** @} */ - -#ifdef __cplusplus -} -#endif/*__cplusplus*/ - - - -/** -@mainpage C port of Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) - -@section intro Introduction - -This library is a C port of the implementation of Limited-memory -Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method written by Jorge Nocedal. -The original FORTRAN source code is available at: -http://www.ece.northwestern.edu/~nocedal/lbfgs.html - -The L-BFGS method solves the unconstrainted minimization problem, - -
-    minimize F(x), x = (x1, x2, ..., xN),
-
- -only if the objective function F(x) and its gradient G(x) are computable. The -Newton's method, which is a well-known algorithm for the optimization, -requires computation or approximation of the inverse of the hessian matrix of -the objective function in order to find the point where the gradient G(X) = 0. -The computational cost for the inverse hessian matrix is expensive especially -when the objective function takes a large number of variables. The L-BFGS -method approximates the inverse hessian matrix efficiently by using information -from last m iterations. This innovation saves the memory storage and -computational time a lot for large-scaled problems. - -Among the various ports of L-BFGS, this library provides several features: -- Optimization with L1-norm (orthant-wise L-BFGS): - In addition to standard minimization problems, the library can minimize - a function F(x) combined with L1-norm |x| of the variables, - {F(x) + C |x|}, where C is a constant scalar parameter. This feature is - useful for estimating parameters of log-linear models with L1-regularization. -- Clean C code: - Unlike C codes generated automatically by f2c (Fortran 77 into C converter), - this port includes changes based on my interpretations, improvements, - optimizations, and clean-ups so that the ported code would be well-suited - for a C code. In addition to comments inherited from the original code, - a number of comments were added through my interpretations. -- Callback interface: - The library receives function and gradient values via a callback interface. - The library also notifies the progress of the optimization by invoking a - callback function. In the original implementation, a user had to set - function and gradient values every time the function returns for obtaining - updated values. -- Thread safe: - The library is thread-safe, which is the secondary gain from the callback - interface. -- Cross platform. The source code can be compiled on Microsoft Visual - Studio 2005, GNU C Compiler (gcc), etc. -- Configurable precision: A user can choose single-precision (float) - or double-precision (double) accuracy by changing ::LBFGS_FLOAT macro. -- SSE/SSE2 optimization: - This library includes SSE/SSE2 optimization (written in compiler intrinsics) - for vector arithmetic operations on Intel/AMD processors. The library uses - SSE for float values and SSE2 for double values. The SSE/SSE2 optimization - routine is disabled by default; compile the library with __SSE__ symbol - defined to activate the optimization routine. - -This library is used by the -CRFsuite project. - -@section download Download - -- Source code - -libLBFGS is distributed under the term of the -MIT license. - -@section changelog History -- Version 1.3 (2007-12-16): - - An API change. An argument was added to lbfgs() function to receive the - final value of the objective function. This argument can be set to - \c NULL if the final value is unnecessary. - - Fixed a null-pointer bug in the sample code (reported by Takashi Imamichi). - - Added build scripts for Microsoft Visual Studio 2005 and GCC. - - Added README file. -- Version 1.2 (2007-12-13): - - Fixed a serious bug in orthant-wise L-BFGS. - An important variable was used without initialization. -- Version 1.1 (2007-12-01): - - Implemented orthant-wise L-BFGS. - - Implemented lbfgs_parameter_init() function. - - Fixed several bugs. - - API documentation. -- Version 1.0 (2007-09-20): - - Initial release. - -@section api Documentation - -- @ref liblbfgs_api "libLBFGS API" - -@section sample Sample code - -@include main.c - -@section ack Acknowledgements - -The L-BFGS algorithm is described in: - - Jorge Nocedal. - Updating Quasi-Newton Matrices with Limited Storage. - Mathematics of Computation, Vol. 35, No. 151, pp. 773--782, 1980. - - Dong C. Liu and Jorge Nocedal. - On the limited memory BFGS method for large scale optimization. - Mathematical Programming B, Vol. 45, No. 3, pp. 503-528, 1989. - -The line search algorithms used in this implementation are described in: - - John E. Dennis and Robert B. Schnabel. - Numerical Methods for Unconstrained Optimization and Nonlinear - Equations, Englewood Cliffs, 1983. - - 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. - -This library also implements Orthant-Wise Limited-memory Quasi-Newton (OW-LQN) -method presented in: - - Galen Andrew and Jianfeng Gao. - Scalable training of L1-regularized log-linear models. - In Proceedings of the 24th International Conference on Machine - Learning (ICML 2007), pp. 33-40, 2007. - -Finally I would like to thank the original author, Jorge Nocedal, who has been -distributing the effieicnt and explanatory implementation in an open source -licence. - -@section reference Reference - -- L-BFGS by Jorge Nocedal. -- OWL-QN by Galen Andrew. -- C port (via f2c) by Taku Kudo. -- C#/C++/Delphi/VisualBasic6 port in ALGLIB. -- Computational Crystallography Toolbox includes - scitbx::lbfgs. -*/ - -#endif/*__LBFGS_H__*/ diff --git a/t/01-parameter.t b/t/01-parameter.t index 6cb3a77..73b58fa 100644 --- a/t/01-parameter.t +++ b/t/01-parameter.t @@ -50,7 +50,7 @@ is_deeply [ 6, 0, - 20 + 40 ], $__;