bfgs算法c语言,机器学习算法实现解析——liblbfgs之L-BFGS算法
在博文“优化算法——拟牛顿法之L-BFGS算法”中,已经对L-BFGS的算法原理做了详细的介绍,本文主要就开源代码liblbfgs重新回顾L-BFGS的算法原理以及具体的实现过程,在L-BFGS算法中包含了处理L1正则的OWL-QN算法,对于OWL-QN算法的详细原理,可以参见博文“优化算法——OWL-QN”。
1、liblbfgs简介
liblbfgs是L-BFGS算法的C语言实现,用于求解非线性优化问题。
下载链接(见上面的主页链接):
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, ¶m);// gp是梯度
} else {// 包含L1正则
ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, ¶m);// 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;
}
收敛的判断方法为:
如果上式满足,则表示已满足收敛条件。
目标函数值是否有足够大的下降(最小问题)
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次的目标函数值。计算的方法为:
是否达到最大的迭代次数// 已达到最大的迭代次数
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矩阵为:
L-BFGS的具体原理可以参见“优化算法——拟牛顿法之L-BFGS算法”。
在上述过程中,第一个循环计算出倒数第m代时的下降方向,第二个阶段利用上面计算出的方法迭代计算出当前的下降方向。
根据上述的流程,开始拟合Hessian矩阵:
// 更新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);// 计算alpha
it->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中涉及到大量的线搜索的策略,线搜索的策略主要作用是找到最优的步长。我们将在下一个主题中进行详细的分析。
补充
1、回调函数
回调函数就是一种利用函数指针进行函数调用的过程。回调函数的好处是具体的计算过程以函数指针的形式传递给其调用处,这样可以较方便地对调用函数进行扩展。
假设有个print_result函数,需要输出两个int型数的和,那么直接写即可,如果需要改成差,那么得重新修改;如果在print_result函数的参数中传入一个函数指针,具体的计算过程在该函数中实现,那么就可以在不改变print_result函数的条件下对功能进行扩充,如下面的例子:
frame.h
#include
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
#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;
}
参考文献
bfgs算法c语言,机器学习算法实现解析——liblbfgs之L-BFGS算法相关推荐
- 狄斯奎诺算法 c语言,图的邻接表实现迪杰斯特拉算法(C语言).doc
图的邻接表实现迪杰斯特拉算法(C语言) /*迪杰斯特拉算法(狄斯奎诺算法)解决的是从源点到其它所有顶点的最短路径问题*/ //算法实现: #include #include #define MAX 2 ...
- 狄斯奎诺算法 c语言,图的邻接表实现迪杰斯特拉算法(C语言)
图的邻接表实现迪杰斯特拉算法(C语言). 迪杰斯特拉算法(狄斯奎诺算法)解决的是从源点到其它所有顶点的最短路径问题. 图的邻接表实现迪杰斯特拉算法(C语言) /*迪杰斯特拉算法(狄斯奎诺算法)解决的是 ...
- 活性边表算法c语言,《计算机图形学》有序边表填充算法.doc
PAGE PAGE 8 实 验 报 告 实验目的 掌握有序边表算法填充多边形区域: 理解多边形填充算法的意义: 增强C语言编程能力. 算法原理介绍 根据多边形内部点的连续性知:一条扫描线与多边形的交点 ...
- sunday算法c语言实现,C / C++学习笔记:实现Sunday算法
Sunday算法 Sunday 算法于 1990 年 Daniel M.Sunday 提出的字符串模式匹配.其效率在匹配随机的字符串时比其他匹配算法还要更快.Sunday 算法的实现可比 KMP,BM ...
- czt算法c语言实现,基--2按频率抽取的FFT算法Decimation-in-Frequency(DIF).ppt
基--2按频率抽取的FFT算法Decimation-in-Frequency(DIF) 第四节基--2按频率抽取的FFT算法Decimation-in-Frequency(DIF)(Sander-Tu ...
- 单像素骨架提取算法c语言实现,【图像】骨架提取与分水岭算法
1.骨架提取 骨架提取,也叫二值图像细化.这种算法能将一个连通区域细化成一个像素的宽度,用于特征提取和目标拓扑表示. morphology子模块提供了两个函数用于骨架提取,分别是Skeletonize ...
- C++算法篇:DFS超详细解析(2)--- tarjan算法求无向图割边
<<<上一篇 系列文章目录 ①:无向图基本概念 ②:tarjan算法求无向图割边 前言 第一次写算法,讲得肯不透彻,有误还请指教awa 文章目录 系列文章目录 一.回顾 二.tarj ...
- mallat算法 c语言,基于STM32F4的小波分解(Mallat算法)程序说明
一.主要思路 原始信号:OrgSig 信号长度:DWT_SIG_LEN 小波分解层数:N 与MATLAB类似,小波分解后产生2个数组DWT_L和DWT_C,但定义与MATLAB不同.定义如下: DWT ...
- aloha算法c语言代码,论文对最简单的反碰撞算法ALOHA算法进行了研究,在识别时间...
论文对最简单的反碰撞算法ALOHA算法进行了研究,在识别时间和重发次数之间作一下折衷,确定如何选择退避时间. 相关句子 1.确定如何选择这三个之一超出了本文的范围. 2.时间有限,精力有限,金钱有限. ...
最新文章
- hive python udf_python udf方法
- 数据结构源码笔记(C语言):英文单词按字典序排序的基数排序
- TortoiseSVN操作
- ubyntu 链接mysql_ubuntu mysql远程连接
- 2020,这些前沿技术成全球关注热点
- for vue 一行2列_前端开发面试问什么?vue面试中经常问到的问题?用vue想拿20k,面试题要这样答!...
- 基于hadoop构建对象存储系统_基于Hadoop企业私有云存储平台的构建
- Java编程:多路查找树
- 阶段3 3.SpringMVC·_03.SpringMVC常用注解_8 SessionAttributes注解
- Android框架揭秘电子书pdf下载
- 项目管理:项目质量管理
- mx播放器有没有投屏功能_这个播放器真是太强大了!
- 【HTML】铺满背景图片
- HDU 6194 后缀数组+单调栈
- 2021.7.15 jzoj题解与反思
- UI设计年薪20W?为什么UI设计能这么火呢?
- 吴恩达Coursera深度学习课程 deeplearning.ai (4-4) 人脸识别和神经风格转换--编程作业
- My personal website:http://47.94.240.229:8080/yjh/project/
- 滚动代码Marquee详解(html滚动显示文字)
- 给你最强的秘籍IE7最新使用技巧荟萃(转)
热门文章
- java -p_说说javap命令
- Python命令行创建工具包——Click(1)
- 文字转换成语音的工具叫什么?这几款软件你确定不试试吗?
- 写作技巧~100段作文排比句(21-40段),考试一定用得上,赶紧收藏!
- ffmpeg播放器实现详解 - 音频播放
- 我一定要赚好多好的钱,孝敬父母
- sriov-network-device-plugin 代码分析
- 树莓派实现云直播系统
- 计算机基础知识用语,电脑基础知识之常见术语
- 网易员工酒吧蹦迪偶遇丁老板,点了一桌价值50w的香槟,愣是一口没喝,原因竟是......