敞开成长之旅！这是我参加「日新方案 2 月更文应战」的第 5 天，点击检查活动概况

在博文“优化算法——拟牛顿法之L-BFGS算法”中，已经对L-BFGS的算法原理做了详细的介绍，本文首要就开源代码liblbfgs从头回顾L-BFGS的算法原理以及详细的完成进程，在L-BFGS算法中包括了处理L1正则的OWL-QN算法，关于OWL-QN算法的详细原理，能够拜见博文“优化算法——OWL-QN”。

1. liblbfgs概述

liblbfgs是根据C语言完成的L-BFGS算法库，用于求解非线性优化问题。能够经过liblbfgs的主页（www.chokkan.org/software/li…）查询到对liblbfgs模块的介绍。其代码能够经过以下的链接下载：

用于Linux渠道 github.com/downloads/c…
用于Windows渠道 github.com/chokkan/lib…

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$

假如上式满意，则表示已满意收敛条件。

方针函数值是否有满足大的下降（最小问题）

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$

是否到达最大的迭代次数

// 已到达最大的迭代次数
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}

其间， $yk=▽f(xk+1)−▽f(xk)y_k=\bigtriangledown f\left ( x_{k+1} \right )-\bigtriangledown f\left ( x_{k} \right )$ ， $s_k=x_{k+1}-x_{k}$ 。在核算的进程中，需求不断的核算和存储前史的Hessian矩阵，在L-BFGS算法，期望只保存最近的 $m$ 次迭代信息，便能够拟合Hessian矩阵。在L-BFGS算法中，不再保存完好的 $H_k$ ，而是存储向量序列 ${sk}\left \{ s_k \right \}$ 和 ${yk}\left \{ y_k \right \}$ ，需求矩阵时 $H_k$ ，运用向量序列 ${sk}\left \{ s_k \right \}$ 和 ${yk}\left \{ y_k \right \}$ 核算就能够得到，而向量序列 ${sk}\left \{ s_k \right \}$ 和 ${yk}\left \{ y_k \right \}$ 也不是一切都要保存，只要保存最新的 $m$ 步向量即可。其详细的核算办法为：

L-BFGS的详细原理能够拜见“优化算法——拟牛顿法之L-BFGS算法”。

在上述进程中，第一个循环核算出倒数第 $m$ 代时的下降方向，第二个阶段运用上面核算出的办法迭代核算出当时的下降方向。

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

核算向量序列 ${sk}\left \{ s_k \right \}$ 和 ${yk}\left \{ y_k \right \}$

// 更新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中涉及到很多的线查找的战略，线查找的战略首要效果是找到最优的步长。咱们将在下一个主题中进行详细的分析。

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

敞开成长之旅！这是我参加「日新方案 2 月更文应战」的第 5 天，点击检查活动概况

liblbfgs中L-BFGS算法的实现

1. liblbfgs概述

2. liblbfgs源码解析

2.1. 源码的首要结构

2.2. 东西`arithmetic_ansi.h`

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

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

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

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

2.3.4. 参数的声明和初始化

2.3.5. 循环求解的进程

2.3.6. 循环的停止条件

2.3.7. 拟合Hessian矩阵

2.3.8. 线查找战略

3. 弥补

3.1. 回调函数

参考文献

作者信息

liblbfgs中L-BFGS算法的实现

1. liblbfgs概述

2. liblbfgs源码解析

2.1. 源码的首要结构

2.2. 东西arithmetic_ansi.h

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

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

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

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

2.3.4. 参数的声明和初始化

2.3.5. 循环求解的进程

2.3.6. 循环的停止条件

2.3.7. 拟合Hessian矩阵

2.3.8. 线查找战略

3. 弥补

3.1. 回调函数

参考文献

相关文章

SkipList 基本原理和Java实现

论文复现丨基于ModelArts实现Text2SQL

Android搞机之编译Android13并在模拟器上运行

Glide 源码阅读笔记（六）

作者信息

2.2. 东西`arithmetic_ansi.h`