在博文“优化算法——拟牛顿法之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;
}

收敛的判断方法为：

$$\frac{\left|g\left(x\right)\right|}{max\left(1,\left|x\right|\right)}<\epsilon$$

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

• 目标函数值是否有足够大的下降（最小问题）

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次的目标函数值。计算的方法为：

$$\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矩阵为：

$$H_{k+1}={\left(I-\frac{s_k{y}_k^T}{y_k^T{s}_k}\right)}^T{H}_k\left(I-\frac{y_k{s}_k^T}{y_k^T{s}_k}\right)+\frac{s_k{s}_k^T}{y_k^T{s}_k}$$

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

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

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

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

• 计算向量序列 { s k } \left \{ s_k \right \} {sk​}和 { y k } \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