在博文“优化算法——拟牛顿法之L-BFGS算法”中,已经对L-BFGS的算法原理做了详细的介绍,本文主要就开源代码liblbfgs重新回顾L-BFGS的算法原理以及具体的实现过程,在L-BFGS算法中包含了处理L1正则的OWL-QN算法,对于OWL-QN算法的详细原理,可以参见博文“优化算法——OWL-QN”。
liblbfgs是基于C语言实现的L-BFGS算法库,用于求解非线性优化问题。可以通过liblbfgs的主页(http://www.chokkan.org/software/liblbfgs/)查询到对liblbfgs模块的介绍。其代码可以通过以下的链接下载:
在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
:测试样例arithmetic_ansi.h
在arithmetic_ansi.h
文件中,主要是对向量(vector)的一些操作,这部分的程序代码比较简单,在这里简单对没个函数的作用进行描述,包括:
void* vecalloc(size_t size)
void vecfree(void *memblock)
void vecset(lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
void veccpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
void vecncpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
void vecadd(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
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)
在liblbfgs中,有很多利用汇编语言优化的代码,这里暂且不考虑这些优化的代码,对于这些优化的代码,作者提供了基本的实现方式。
函数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);
}
函数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参数的结构体,其具体的代码如下所示:
作者在这部分代码中的注释写得特别详细,从这些注释中可以学习到很多调参时的重要信息。
函数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
中,其基本的过程为:
在初始化阶段涉及到很多参数的声明,接下来将详细介绍每一个参数的作用。
循环的求解过程从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中的梯度。
在选择了最优步长过程中,会同时对变量进行更新,第二步即是判断此时的更新是否满足终止条件,终止条件分为以下三类:
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矩阵。
在BFGS算法(优化算法——拟牛顿法之BFGS算法)中,其Hessian矩阵为:
其中,
,
。在计算的过程中,需要不断的计算和存储历史的Hessian矩阵,在L-BFGS算法,希望只保留最近的
次迭代信息,便能够拟合Hessian矩阵。在L-BFGS算法中,不再保存完整的
,而是存储向量序列
和
,需要矩阵时
,使用向量序列
和
计算就可以得到,而向量序列
和
也不是所有都要保存,只要保存最新的
步向量即可。其具体的计算方法为:
L-BFGS的具体原理可以参见“优化算法——拟牛顿法之L-BFGS算法”。
在上述过程中,第一个循环计算出倒数第
代时的下降方向,第二个阶段利用上面计算出的方法迭代计算出当前的下降方向。
根据上述的流程,开始拟合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;
在liblbfgs中涉及到大量的线搜索的策略,线搜索的策略主要作用是找到最优的步长。我们将在下一个主题中进行详细的分析。
回调函数就是一种利用函数指针进行函数调用的过程。回调函数的好处是具体的计算过程以函数指针的形式传递给其调用处,这样可以较方便地对调用函数进行扩展。
假设有个print_result
函数,需要输出两个int型数的和,那么直接写即可,如果需要改成差,那么得重新修改;如果在print_result
函数的参数中传入一个函数指针,具体的计算过程在该函数中实现,那么就可以在不改变print_result
函数的条件下对功能进行扩充,如下面的例子:
#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));
}
#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;
}