在博文“优化算法——拟牛顿法之L-BFGS算法”中，已经对L-BFGS的算法原理做了详细的介绍，本文主要就开源代码liblbfgs重新回顾L-BFGS的算法原理以及具体的实现过程，在L-BFGS算法中包含了处理L1正则的OWL-QN算法，对于OWL-QN算法的详细原理，可以参见博文“优化算法——OWL-QN”。

1. liblbfgs概述

liblbfgs是基于C语言实现的L-BFGS算法库，用于求解非线性优化问题。可以通过liblbfgs的主页（http://www.chokkan.org/software/liblbfgs/）查询到对liblbfgs模块的介绍。其代码可以通过以下的链接下载：

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

2. liblbfgs源码解析

2.1. 源码的主要结构

在liblbfgs中，主要的代码包括

  • liblbfgs-1.10/include/lbfgs.h：头文件
  • liblbfgs-1.10/lib/lbfgs.c：具体的实现
  • liblbfgs-1.10/lib/arithmetic_ansi.h(另两个arithmetic_sse_double.h和arithmetic_sse_float.h是两个汇编编写的等价形式)：相当于一些工具
  • liblbfgs-1.10/sample/sample.c：测试样例

2.2. 工具arithmetic_ansi.h

在arithmetic_ansi.h文件中，主要是对向量（vector）的一些操作，这部分的程序代码比较简单，在这里简单对没个函数的作用进行描述，包括：

  • 申请size大小的空间，同时对其进行初始化
void* vecalloc(size_t size)
  • 释放空间
void vecfree(void *memblock)
  • 将向量x中的值设置为c
void vecset(lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
  • 将向量x中的值复制到向量y中
void veccpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
  • 取向量x中的每个值的负数，将其放到向量y中
void vecncpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
  • 对向量y中的每个元素增加向量x中对应元素的c倍
void vecadd(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
  • 计算向量x和向量y的差
void vecdiff(lbfgsfloatval_t *z, const lbfgsfloatval_t *x, const lbfgsfloatval_t *y, const int n)
  • 向量与常数的积
void vecscale(lbfgsfloatval_t *y, const lbfgsfloatval_t c, const int n)
  • 向量的外积
void vecmul(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
  • 向量的点积
void vecdot(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const lbfgsfloatval_t *y, const int n)
  • 向量的点积的开方
void vec2norm(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n)
  • 向量的点积的开方的倒数
void vec2norminv(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n)

2.3. L-BFGS算法的主要函数

在liblbfgs中，有很多利用汇编语言优化的代码，这里暂且不考虑这些优化的代码，对于这些优化的代码，作者提供了基本的实现方式。

2.3.1. 为变量分配和回收内存空间

函数lbfgs_malloc是为优化函数中的变量分配内存空间的函数，其具体形式为：

// 为变量分配空间
lbfgsfloatval_t* lbfgs_malloc(int n)
{// 涉及到汇编的一些知识，暂且不考虑
#if     defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__))n = round_out_variables(n);
#endif/*defined(USE_SSE)*/// 分配n个大小的内存空间return (lbfgsfloatval_t*)vecalloc(sizeof(lbfgsfloatval_t) * n);
}

函数lbfgs_malloc用于为变量分配长度为n的内存空间，其中，lbfgsfloatval_t为定义的变量的类型，其定义为float或者是double类型：

#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

与内存分配对应的是内存的回收，其具体形式为：

// 释放变量的空间
void lbfgs_free(lbfgsfloatval_t *x)
{vecfree(x);
}

2.3.2. L-BFGS中参数的初始化

函数lbfgs_parameter_init提供了L-BFGS默认参数的初始化方法。

其实在L-BFGS的算法过程中也会提供默认的参数的方法，所以该方法有点多余。

// 默认初始化lbfgs的参数
void lbfgs_parameter_init(lbfgs_parameter_t *param)
{memcpy(param, &_defparam, sizeof(*param));// 使用默认的参数初始化
}

函数lbfgs_parameter_init将默认参数_defparam复制到参数param中，lbfgs_parameter_t是L-BFGS参数的结构体，其具体的代码如下所示：

作者在这部分代码中的注释写得特别详细，从这些注释中可以学习到很多调参时的重要信息。

2.3.3. L-BFGS算法的核心过程概述

函数lbfgs是优化算法的核心过程，其函数的参数为：

int n,// 变量的个数
lbfgsfloatval_t *x,// 变量
lbfgsfloatval_t *ptr_fx,// 目标函数值
lbfgs_evaluate_t proc_evaluate,// 计算目标函数值和梯度的回调函数
lbfgs_progress_t proc_progress,// 处理计算过程的回调函数
void *instance,// 数据
lbfgs_parameter_t *_param// L-BFGS的参数

该函数通过调用两个函数proc_evaluate和proc_progress用于计算具体的函数以及处理计算过程中的一些工作。

在函数lbfgs中，其基本的过程为：

2.3.4. 参数的声明和初始化

在初始化阶段涉及到很多参数的声明，接下来将详细介绍每一个参数的作用。

2.3.5. 循环求解的过程

循环的求解过程从for循环开始，在for循环中，首先需要利用线搜索策略进行最优的步长选择，其具体的过程如下所示：

/* Store the current position and gradient vectors. */
// 存储当前的变量值和梯度值
veccpy(xp, x, n);// 将当前的变量值复制到向量xp中
veccpy(gp, g, n);// 将当前的梯度值复制到向量gp中/* Search for an optimal step. */
// 线搜索策略，搜索最优步长
if (param.orthantwise_c == 0.) {// 无L1正则ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, &param);// gp是梯度
} else {// 包含L1正则ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, &param);// pg是伪梯度// 计算伪梯度owlqn_pseudo_gradient(pg, x, g, n,param.orthantwise_c, param.orthantwise_start, param.orthantwise_end);
}
if (ls < 0) {// 已达到终止条件// 由于在线搜索的过程中更新了向量x和向量g，因此从xp和gp中恢复变量值和梯度值/* Revert to the previous point. */veccpy(x, xp, n);veccpy(g, gp, n);ret = ls;goto lbfgs_exit;// 释放资源
}

由于在线搜索的过程中会对变量x以及梯度的向量g修改，因此在进行线搜索之前，先将变量x以及梯度g保存到另外的向量中。参数param.orthantwise_c表示的是L1正则的参数，若为0则不使用L1正则，即使用L-BFGS算法；若不为0，则使用L1正则，即使用OWL-QN算法。

关于线搜索的具体方法，可以参见2.3.6。
对于owlqn_pseudo_gradient函数，可以参见2.3.4

在OWL-QN中，由于在某些点处不存在导数，因此使用伪梯度代替L-BFGS中的梯度。

2.3.6. 循环的终止条件

在选择了最优步长过程中，会同时对变量进行更新，第二步即是判断此时的更新是否满足终止条件，终止条件分为以下三类：

  • 是否收敛
vec2norm(&xnorm, x, n);// 平方和的开方
if (param.orthantwise_c == 0.) {// 非L1正则vec2norm(&gnorm, g, n);
} else {// L1正则vec2norm(&gnorm, pg, n);
}

// 判断是否收敛
if (xnorm < 1.0) xnorm = 1.0;
if (gnorm / xnorm <= param.epsilon) {/* Convergence. */ret = LBFGS_SUCCESS;break;
}

收敛的判断方法为：

∣g(x)∣max(1,∣x∣)<ε\frac{\left | g\left ( x \right ) \right |}{max\left ( 1,\left | x \right | \right )}<\varepsilon max(1,∣x∣)∣g(x)∣​<ε

如果上式满足，则表示已满足收敛条件。

  • 目标函数值是否有足够大的下降（最小问题）
if (pf != NULL) {// 终止条件/* We don't test the stopping criterion while k < past. */// k为迭代次数，只考虑past>k的情况，past是指只保留past次的值if (param.past <= k) {/* Compute the relative improvement from the past. */// 计算函数减小的比例rate = (pf[k % param.past] - fx) / fx;/* The stopping criterion. */// 下降比例是否满足条件if (rate < param.delta) {ret = LBFGS_STOP;break;}}/* Store the current value of the objective function. */// 更新新的目标函数值pf[k % param.past] = fx;
}

在pf中，保存了param.past次的目标函数值。计算的方法为：

fo−fnfn<Δ\frac{f_o-f_n}{f_n}<\Delta fn​fo​−fn​​<Δ

  • 是否达到最大的迭代次数
// 已达到最大的迭代次数
if (param.max_iterations != 0 && param.max_iterations < k+1) {/* Maximum number of iterations. */ret = LBFGSERR_MAXIMUMITERATION;break;
}

如果没有满足终止的条件，那么需要拟合新的Hessian矩阵。

2.3.7. 拟合Hessian矩阵

在BFGS算法（优化算法——拟牛顿法之BFGS算法）中，其Hessian矩阵为：

Hk+1=(I−skykTykTsk)THk(I−ykskTykTsk)+skskTykTskH_{k+1}=\left ( I-\frac{s_ky_k^T}{y_k^Ts_k} \right )^TH_k\left ( I-\frac{y_ks_k^T}{y_k^Ts_k} \right )+\frac{s_ks_k^T}{y_k^Ts_k}Hk+1​=(I−ykT​sk​sk​ykT​​)THk​(I−ykT​sk​yk​skT​​)+ykT​sk​sk​skT​​

其中，yk=▽f(xk+1)−▽f(xk)y_k=\bigtriangledown f\left ( x_{k+1} \right )-\bigtriangledown f\left ( x_{k} \right )yk​=▽f(xk+1​)−▽f(xk​)，sk=xk+1−xks_k=x_{k+1}-x_{k}sk​=xk+1​−xk​。在计算的过程中，需要不断的计算和存储历史的Hessian矩阵，在L-BFGS算法，希望只保留最近的mmm次迭代信息，便能够拟合Hessian矩阵。在L-BFGS算法中，不再保存完整的HkH_kHk​，而是存储向量序列{sk}\left \{ s_k \right \}{sk​}和{yk}\left \{ y_k \right \}{yk​}，需要矩阵时HkH_kHk​，使用向量序列{sk}\left \{ s_k \right \}{sk​}和{yk}\left \{ y_k \right \}{yk​}计算就可以得到，而向量序列{sk}\left \{ s_k \right \}{sk​}和{yk}\left \{ y_k \right \}{yk​}也不是所有都要保存，只要保存最新的mmm步向量即可。其具体的计算方法为：

L-BFGS的具体原理可以参见“优化算法——拟牛顿法之L-BFGS算法”。

在上述过程中，第一个循环计算出倒数第mmm代时的下降方向，第二个阶段利用上面计算出的方法迭代计算出当前的下降方向。

根据上述的流程，开始拟合Hessian矩阵：

  • 计算向量序列{sk}\left \{ s_k \right \}{sk​}和{yk}\left \{ y_k \right \}{yk​}
// 更新s向量和y向量
it = &lm[end];// 初始时，end为0
vecdiff(it->s, x, xp, n);// x - xp，xp为上一代的值，x为当前的值
vecdiff(it->y, g, gp, n);// g - gp，gp为上一代的值，g为当前的值
  • 两个循环
// 通过拟牛顿法计算Hessian矩阵
// L-BFGS的两个循环
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);// 计算alphait->alpha /= it->ys;// 乘以rho/* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */vecadd(d, it->y, -it->alpha, n);// 重新修正d
}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;// 乘以rho/* \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; */
}

其中，ys和yy的计算方法如下所示：

vecdot(&ys, it->y, it->s, n);// 计算点积
vecdot(&yy, it->y, it->y, n);
it->ys = ys;

bound和end的计算方法如下所示：

bound = (m <= k) ? m : k;// 判断是否有足够的m代
++k;
end = (end + 1) % m;

2.3.8. 线搜索策略

在liblbfgs中涉及到大量的线搜索的策略，线搜索的策略主要作用是找到最优的步长。我们将在下一个主题中进行详细的分析。

3. 补充

3.1. 回调函数

回调函数就是一种利用函数指针进行函数调用的过程。回调函数的好处是具体的计算过程以函数指针的形式传递给其调用处，这样可以较方便地对调用函数进行扩展。

假设有个print_result函数，需要输出两个int型数的和，那么直接写即可，如果需要改成差，那么得重新修改；如果在print_result函数的参数中传入一个函数指针，具体的计算过程在该函数中实现，那么就可以在不改变print_result函数的条件下对功能进行扩充，如下面的例子：

  • frame.h
#include <stdio.h>void print_result(int (*get_result)(int, int), int a, int b){printf("the final result is : %d\n", get_result(a, b));
}
  • process.cc
#include <stdio.h>
#include "frame.h"int add(int a, int b){return a + b;
}int sub(int a, int b){return a - b;
}int mul(int a, int b){return a * b;
}int max(int a, int b){return (a>b?a:b);
}int main(){int a = 1;int b = 2;print_result(add, a, b);print_result(sub, a, b);print_result(mul, a, b);print_result(max, a, b);return 1;
}

参考文献

  • 优化算法——拟牛顿法之L-BFGS算法
  • 优化算法——OWL-QN

liblbfgs中L-BFGS算法的实现相关推荐

  1. 讨论帖：比特币中的SHA256算法的实现与标准的SHA256算法实现是否相同？

    近日阅读了比特币源码中与哈希相关的部分,对于其中一些细节还是有不清晰的地方. 于是我写了一个小的测试demo:sha256_test,代码下载 分别测试了三个版本对于SHA-256算法的实现: Bit ...

  2. ROS中，DWA算法的实现

    在ROS中,DWA算法的实现主要涉及到以下几个方面: 机器人运动学模型:DWA算法需要机器人的运动学模型,ROS中提供了很多机器人模型,可以根据实际情况进行选择. 环境地图:DWA算法需要环境地图,R ...

  3. linkedhashmap中关于LRU算法的实现

    //LinkedHashMap的一个构造函数,当参数accessOrder为true时,即会按照访问顺序排序,最近访问的放在最前,最早访问的放在后面 public LinkedHashMap(int ...

  4. PCA（主成分分析法）的理解笔记及算法的实现

    前几天搞定了Open3d库问题后,准备手撕PCA算法突然人麻了.我坚信学习是不断重复的过程,特此做个笔记,欢迎大家评论和交流! 感谢大佬的文章: 1.主成分分析(PCA)原理详解_Microstron ...

  5. 53.垃圾回收算法的实现原理、启动Java垃圾回收、Java垃圾回收过程、垃圾回收中实例的终结、对象什么时候符合垃圾回收的条件、GC Scope 示例程序、GC OutOfMemoryError的示例

    53.垃圾回收算法的实现原理 53.1.目录 53.2.启动Java垃圾回收 53.3.Java垃圾回收过程 53.4.垃圾回收中实例的终结 53.5.对象什么时候符合垃圾回收的条件? 53.5.1. ...

  6. 计算机LCG/PCG/MWC/XorShift等PRNG算法，以及V8中Math.random()、webkit中crypto等随机算法的实现

    计算机LCG/PCG/MWC/XorShift等PRNG算法,以及V8中Math.random().webkit中crypto等随机算法的实现 本文篇幅较长,如想直接看 js 的随机数实现可定位本文E ...

  7. OpenCV中图像旋转(warpAffine)算法的实现过程

    在OpenCV中,目前并没有现成的函数直接用来实现图像旋转,它是用仿射变换函数cv::warpAffine来实现的,此函数目前支持4种插值算法,最近邻.双线性.双三次.兰索斯插值,如果传进去的参数为基 ...

  8. JAVA实现中点画线_实验1-中点画线和Bresenham画线算法的实现

    <实验1-中点画线和Bresenham画线算法的实现>由会员分享,可在线阅读,更多相关<实验1-中点画线和Bresenham画线算法的实现(9页珍藏版)>请在人人文库网上搜索. ...

  9. 游戏中常用的寻路算法的分享（3）：A*算法的实现

    概述 剥除代码,A* 算法非常简单.算法维护两个集合:OPEN 集和 CLOSED 集.OPEN 集包含待检测节点.初始状态,OPEN集仅包含一个元素:开始位置.CLOSED集包含已检测节点.初始状态 ...

最新文章

  1. WCF服务重构实录（上）
  2. python代码规范化_数据标准化方法及其Python代码实现
  3. 关于 printk() 对 spi slave 内核驱动程序的性能影响
  4. 软件架构设计最佳实践
  5. 【IT笔试面试题整理】判断一个树是否是另一个的子树
  6. 使用 .NET Core模板引擎创建自定义的模板和项目
  7. 误操作删除数据文件恢复案例讨论
  8. node连接--MySQL
  9. 正负样本不平衡处理方法总结【转】
  10. Hadoop-HDFS原理及操作（小实验）
  11. SpringBoot+websocket 实现web聊天功能（单聊、保存消息）
  12. 数据分析在银行业应用之欺诈检测
  13. 51单片机实现拼音输入法
  14. 走进MyBatis源码一探Spring扩展点「知识点多多」「扩展点实战系列」- 第449篇
  15. Spring框架的七大模块
  16. 苹果电池测试软件原理,苹果手机电池检测软件有哪些?
  17. DNA-蛋白翻译过程的Python实现
  18. 【SCAU 新生赛】18247 aler的旅游计划 并查集模板题
  19. Lora、zigbee比较
  20. H5头像完整制作，可拖拽缩放，可添加装饰图标（装饰图标支持缩放、旋转、拖拽）

热门文章

  1. 低压铸造水模拟计算机控制系统,低压铸造中充型过程水模拟的研究.doc
  2. 新思路｜保姆式智能化安防方案，开启全新商业管理模式
  3. 【训练题23：中国剩余定理】猜数字 | P3868 [TJOI2009]
  4. charles修改下行参数
  5. 年味贺岁！雅迪“黑科技”为非遗文化注入年轻活力
  6. 什么浏览器最好用,浏览器大排行！
  7. vue实现省市区三级联动地址选择
  8. “换脸术”即将带来的“福利”与恐慌
  9. 360 2018年春招编程题第三题
  10. 安川机器人的变量类型