非线性优化库liblbfgs初探

liblbfgs简介

liblbfgs是L-BFGS算法的C语言实现，用于求解非线性优化问题。

liblbfgs的主页：http://www.chokkan.org/software/liblbfgs/

下载链接（见上面的主页链接）：

• https://github.com/downloads/chokkan/liblbfgs/liblbfgs-1.10.tar.gz 用于Linux平台
• https://github.com/chokkan/liblbfgs 
  用于Windows平台

Linux下用gcc 编译

解压文件，进入目录，可参照 软件包里面的说明文档，命令依次如下：

./configure

make

make install

编译后的目录文件：

C源码文件在lib,include和 sample 三个文件夹里，而主文件在’sample’ 里，其目录如下：

其中没有文件后缀名的文件sample就是最终的可执行文件，运行： 

Windows下用visual studio编译

只需要把几个C/C++源码文件加入到项目即可，其余的文件一律弃之不用，如下图所示的VS项目文件的目录，我们只要 4 个头文件 arithmetic_ansi.h， arithmetic_sse_double.h， arithmetic_sse_float.h，lbfgs.h 和三个源文件 
lbfgs.c， sample.c， sample.cpp。

主函数文件（包含main函数）是 sample.c（用C语言写的），和 sample.cpp（用C++写成，封装成一个类）。一个VS项目里只能有一个main函数，所以每次运行要禁用二者中的一个。

需要注意的是头文件中的 arithmetic_ansi.h， arithmetic_sse_double.h， arithmetic_sse_float.h，三个头文件只能用一个，选择哪个由’lbfgs.cpp’中一段代码控制，该段代码如下：

默认情况下 USE_SSE、__SSE2__、__SSE__ 都没有定义，L-BFGS_FLOAT 在 lbfgs.h中定义。由此可知默认情况下采用的是 arithmetic_ansi.h 。

若要采用 arithmetic_sse_double.h，经查L-BFGS_FLOAT 被定义为 64，所以再定义其它两个变量，如下：

SSE2到底是什么东西呢？插入一段相关介绍：

SSE2(Streaming SIMD Extensions 2)：指令集，扩展了SSE指令集，并可完全取代MMX。SSE2指令集是Intel公司在SSE指令集的基础上发展起来的。相比于SSE，SSE2使用了144个新增指令，扩展了MMX技术和SSE技术，提高了诸如MPEG-2、MP3、3D图形等应用程序的运行性能。在整数处理方面，随MMX技术引进的SIMD整数指令从64位扩展到了128 位，使SIMD整数类型操作的执行效率成倍提高；在浮点数处理方面，双精度浮点SIMD指令允许以 SIMD格式同时执行两个浮点操作，提供双精度操作支持有助于加速内容创建、财务、工程和科学应用。除SSE2指令之外，最初的SSE指令也得到增强，通过支持多种数据类型(例如双字、四字)的算术运算，支持灵活、动态范围更广的计算功能。

可见引入 arithmetic_sse_double.h 或 arithmetic_sse_float.h 是为了提高浮点数的运算速度，改善程序性能。

main函数分析

main函数里给出了一个用L-BFGS算法进行非线性优化的例子，贴出main函数的代码：

<code class="language-c hljs  has-numbering" style="display: block; padding: 0px; color: inherit; box-sizing: border-box; font-family: 'Source Code Pro', monospace;font-size:undefined; white-space: pre; border-top-left-radius: 0px; border-top-right-radius: 0px; border-bottom-right-radius: 0px; border-bottom-left-radius: 0px; word-wrap: normal; background: transparent;"><span class="hljs-preprocessor" style="color: rgb(68, 68, 68); box-sizing: border-box;">#include <stdio.h></span><span class="hljs-preprocessor" style="color: rgb(68, 68, 68); box-sizing: border-box;">#include "lbfgs.h"</span><span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">//注意lbfgsfloatval_t在lbfgs.h里定义，32位机上是float，64位机上是double。</span><span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">//evaluate定义要优化的非线性函数fx及其梯度g</span><span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">static</span> lbfgsfloatval_t evaluate(     <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">void</span> *instance,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t *x,    lbfgsfloatval_t *g,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> n,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t step    ){    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> i;    lbfgsfloatval_t fx = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0.0</span>;    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> (i = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;i < n;i += <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>) {        lbfgsfloatval_t t1 = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0</span> - x[i];        lbfgsfloatval_t t2 = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">10.0</span> * (x[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] - x[i] * x[i]);        g[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">20.0</span> * t2;        g[i] = -<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2.0</span> * (x[i] * g[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] + t1);        fx += t1 * t1 + t2 * t2;    }    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">return</span> fx;}<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">//progress函数每次迭代结束就执行一次，打印迭代信息</span><span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">static</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> progress(    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">void</span> *instance,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t *x,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t *g,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t fx,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t xnorm,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t gnorm,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">const</span> lbfgsfloatval_t step,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> n,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> k,    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> ls    ){    <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"Iteration %d:\n"</span>, k);    <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"  fx = %f,"</span>, fx);    <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"  xnorm = %f, gnorm = %f, step = %f\n"</span>, xnorm, gnorm, step);    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">return</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;}<span class="hljs-preprocessor" style="color: rgb(68, 68, 68); box-sizing: border-box;">#define N   100 <span class="hljs-comment" style="box-sizing: border-box;">//参数个数</span></span><span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> main(<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> argc, <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">char</span> *argv[]){    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">int</span> i, ret = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;    lbfgsfloatval_t fx;    lbfgsfloatval_t *x = lbfgs_malloc(N); <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">//x是长度为N的数组</span>    lbfgs_parameter_t param;    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">if</span> (x == NULL) {        <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"ERROR: Failed to allocate a memory block for variables.\n"</span>);        <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">return</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>;    }    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">/* Initialize the variables. */</span>    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> (i = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;i < N;i += <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>) {        x[i] = -<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.2</span>;        x[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0</span>;    }    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">/* Initialize the parameters for the L-BFGS optimization. */</span>    lbfgs_parameter_init(&param);    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">/*param.linesearch = LBFGS_LINESEARCH_BACKTRACKING;*/</span>    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">/*        Start the L-BFGS optimization; this will invoke the callback functions        evaluate() and progress() when necessary.     */</span>    ret = lbfgs(N, x, &fx, evaluate, progress, NULL, &param);    <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"算法收敛后的函数取值  fx = %f\n"</span>, fx, x[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>], x[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>]);    <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"最终的解（长度为100个浮点数的数组）："</span>);    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span>(i=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;i<N;i++) {        <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">if</span> (i%<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">10</span> == <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>) {            <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"\n"</span>);        }        <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">printf</span>(<span class="hljs-string" style="color: rgb(0, 136, 0); box-sizing: border-box;">"%2.2lf "</span>,x[i]);    }    lbfgs_free(x); <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">//释放参数数组所占的内存</span>    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">return</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;}</code><ul class="pre-numbering" style="box-sizing: border-box; position: absolute; width: 50px; top: 0px; left: 0px; margin: 0px; padding: 6px 0px 40px; border-right-width: 1px; border-right-style: solid; border-right-color: rgb(221, 221, 221); list-style: none; text-align: right; background-color: rgb(238, 238, 238);"><li style="box-sizing: border-box; padding: 0px 5px;">1</li><li style="box-sizing: border-box; padding: 0px 5px;">2</li><li style="box-sizing: border-box; padding: 0px 5px;">3</li><li style="box-sizing: border-box; padding: 0px 5px;">4</li><li style="box-sizing: border-box; padding: 0px 5px;">5</li><li style="box-sizing: border-box; padding: 0px 5px;">6</li><li style="box-sizing: border-box; padding: 0px 5px;">7</li><li style="box-sizing: border-box; padding: 0px 5px;">8</li><li style="box-sizing: border-box; padding: 0px 5px;">9</li><li style="box-sizing: border-box; padding: 0px 5px;">10</li><li style="box-sizing: border-box; padding: 0px 5px;">11</li><li style="box-sizing: border-box; padding: 0px 5px;">12</li><li style="box-sizing: border-box; padding: 0px 5px;">13</li><li style="box-sizing: border-box; padding: 0px 5px;">14</li><li style="box-sizing: border-box; padding: 0px 5px;">15</li><li style="box-sizing: border-box; padding: 0px 5px;">16</li><li style="box-sizing: border-box; padding: 0px 5px;">17</li><li style="box-sizing: border-box; padding: 0px 5px;">18</li><li style="box-sizing: border-box; padding: 0px 5px;">19</li><li style="box-sizing: border-box; padding: 0px 5px;">20</li><li style="box-sizing: border-box; padding: 0px 5px;">21</li><li style="box-sizing: border-box; padding: 0px 5px;">22</li><li style="box-sizing: border-box; padding: 0px 5px;">23</li><li style="box-sizing: border-box; padding: 0px 5px;">24</li><li style="box-sizing: border-box; padding: 0px 5px;">25</li><li style="box-sizing: border-box; padding: 0px 5px;">26</li><li style="box-sizing: border-box; padding: 0px 5px;">27</li><li style="box-sizing: border-box; padding: 0px 5px;">28</li><li style="box-sizing: border-box; padding: 0px 5px;">29</li><li style="box-sizing: border-box; padding: 0px 5px;">30</li><li style="box-sizing: border-box; padding: 0px 5px;">31</li><li style="box-sizing: border-box; padding: 0px 5px;">32</li><li style="box-sizing: border-box; padding: 0px 5px;">33</li><li style="box-sizing: border-box; padding: 0px 5px;">34</li><li style="box-sizing: border-box; padding: 0px 5px;">35</li><li style="box-sizing: border-box; padding: 0px 5px;">36</li><li style="box-sizing: border-box; padding: 0px 5px;">37</li><li style="box-sizing: border-box; padding: 0px 5px;">38</li><li style="box-sizing: border-box; padding: 0px 5px;">39</li><li style="box-sizing: border-box; padding: 0px 5px;">40</li><li style="box-sizing: border-box; padding: 0px 5px;">41</li><li style="box-sizing: border-box; padding: 0px 5px;">42</li><li style="box-sizing: border-box; padding: 0px 5px;">43</li><li style="box-sizing: border-box; padding: 0px 5px;">44</li><li style="box-sizing: border-box; padding: 0px 5px;">45</li><li style="box-sizing: border-box; padding: 0px 5px;">46</li><li style="box-sizing: border-box; padding: 0px 5px;">47</li><li style="box-sizing: border-box; padding: 0px 5px;">48</li><li style="box-sizing: border-box; padding: 0px 5px;">49</li><li style="box-sizing: border-box; padding: 0px 5px;">50</li><li style="box-sizing: border-box; padding: 0px 5px;">51</li><li style="box-sizing: border-box; padding: 0px 5px;">52</li><li style="box-sizing: border-box; padding: 0px 5px;">53</li><li style="box-sizing: border-box; padding: 0px 5px;">54</li><li style="box-sizing: border-box; padding: 0px 5px;">55</li><li style="box-sizing: border-box; padding: 0px 5px;">56</li><li style="box-sizing: border-box; padding: 0px 5px;">57</li><li style="box-sizing: border-box; padding: 0px 5px;">58</li><li style="box-sizing: border-box; padding: 0px 5px;">59</li><li style="box-sizing: border-box; padding: 0px 5px;">60</li><li style="box-sizing: border-box; padding: 0px 5px;">61</li><li style="box-sizing: border-box; padding: 0px 5px;">62</li><li style="box-sizing: border-box; padding: 0px 5px;">63</li><li style="box-sizing: border-box; padding: 0px 5px;">64</li><li style="box-sizing: border-box; padding: 0px 5px;">65</li><li style="box-sizing: border-box; padding: 0px 5px;">66</li><li style="box-sizing: border-box; padding: 0px 5px;">67</li><li style="box-sizing: border-box; padding: 0px 5px;">68</li><li style="box-sizing: border-box; padding: 0px 5px;">69</li><li style="box-sizing: border-box; padding: 0px 5px;">70</li><li style="box-sizing: border-box; padding: 0px 5px;">71</li><li style="box-sizing: border-box; padding: 0px 5px;">72</li><li style="box-sizing: border-box; padding: 0px 5px;">73</li><li style="box-sizing: border-box; padding: 0px 5px;">74</li><li style="box-sizing: border-box; padding: 0px 5px;">75</li><li style="box-sizing: border-box; padding: 0px 5px;">76</li><li style="box-sizing: border-box; padding: 0px 5px;">77</li><li style="box-sizing: border-box; padding: 0px 5px;">78</li><li style="box-sizing: border-box; padding: 0px 5px;">79</li><li style="box-sizing: border-box; padding: 0px 5px;">80</li><li style="box-sizing: border-box; padding: 0px 5px;">81</li><li style="box-sizing: border-box; padding: 0px 5px;">82</li><li style="box-sizing: border-box; padding: 0px 5px;">83</li><li style="box-sizing: border-box; padding: 0px 5px;">84</li><li style="box-sizing: border-box; padding: 0px 5px;">85</li><li style="box-sizing: border-box; padding: 0px 5px;">86</li><li style="box-sizing: border-box; padding: 0px 5px;">87</li><li style="box-sizing: border-box; padding: 0px 5px;">88</li><li style="box-sizing: border-box; padding: 0px 5px;">89</li><li style="box-sizing: border-box; padding: 0px 5px;">90</li><li style="box-sizing: border-box; padding: 0px 5px;">91</li></ul>

非线性函数f(x) 及其 梯度 g(x) 的定义如下：

<code class="language-c hljs  has-numbering" style="display: block; padding: 0px; color: inherit; box-sizing: border-box; font-family: 'Source Code Pro', monospace;font-size:undefined; white-space: pre; border-top-left-radius: 0px; border-top-right-radius: 0px; border-bottom-right-radius: 0px; border-bottom-left-radius: 0px; word-wrap: normal; background: transparent;"><span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> (i = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;i < n;i += <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>) {    lbfgsfloatval_t t1 = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0</span> - x[i];    lbfgsfloatval_t t2 = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">10.0</span> * (x[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] - x[i] * x[i]);    g[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">20.0</span> * t2;    g[i] = -<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2.0</span> * (x[i] * g[i+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] + t1);    fx += t1 * t1 + t2 * t2;}</code><ul class="pre-numbering" style="box-sizing: border-box; position: absolute; width: 50px; top: 0px; left: 0px; margin: 0px; padding: 6px 0px 40px; border-right-width: 1px; border-right-style: solid; border-right-color: rgb(221, 221, 221); list-style: none; text-align: right; background-color: rgb(238, 238, 238);"><li style="box-sizing: border-box; padding: 0px 5px;">1</li><li style="box-sizing: border-box; padding: 0px 5px;">2</li><li style="box-sizing: border-box; padding: 0px 5px;">3</li><li style="box-sizing: border-box; padding: 0px 5px;">4</li><li style="box-sizing: border-box; padding: 0px 5px;">5</li><li style="box-sizing: border-box; padding: 0px 5px;">6</li><li style="box-sizing: border-box; padding: 0px 5px;">7</li></ul>

可知输入参数 有100个，即数组 x，函数f(x)定义如下

f(x)=∑i=0,2,4,...,N−2(1−xi)2+(10(xi+1−x2i))2

其奇数位置(i=1,3,5,...)的参数的梯度 

200(xi+1−x2i)

其偶数位置 (i=0,2,4,6,...)的参数的梯度 

−2(1−xi)−400∗xi∗(xi+1−x2i)

python 封装liblbfgs

• pylbfgs 用Python封装的 liblbfgs 
  https://github.com/larsmans/pylbfgs
• PyLBFGS 
  https://github.com/sseemayer/PyLBFGS/blob/master/test.py