L-BFGS 解优化问题

针对无约束优化问题

$\min_x f(x),$

L-BFGS 在拟牛顿法 BFGS 迭代格式的基础上进行修改，用以解决大规模问题的存储和计算困难。对于拟牛顿法中的迭代方向 $d^k=H^k\nabla f(x^k)$ 考虑利用递归展开的方式进行求解。

首先，对于 BFGS 迭代格式， $H^{k+1}=(V^k)^\top H^kV^k+\rho_ks^k(s^k)^\top$ ，其中 $\displaystyle\rho_k = \frac{1}{(y^k)^\top s^k},$ $V^k=I-\rho_k y^k(s^k)^\top$ 。将其递归地展开得到

$\displaystyle\begin{array}{rl} -H^k\nabla f(x^k)= & -(V^{k-m}\cdots V^{k-1})^\top H^{k-m}(V^{k-m}\cdots V^{k-1})\nabla f(x^k)\\ & -\rho_{k-m}(V^{k-m+1}\cdots V^{k-1})^\top s^{k-m}(s^{k-m})^\top (V^{k-m+1}\cdots V^{k-1})\nabla f(x^k)\\ & -\rho_{k-m+1}(V^{k-m+2}\cdots V^{k-1})^\top s^{k-m+1} (s^{k-m+1})^\top(V^{k-m+2}\cdots V^{k-1})\nabla f(x^k)\\ &-\cdots\\ &-\rho_{k-1}s^{k-1}(s^{k-1^\top})\nabla f(x^k).\end{array}$

我们只需对其中的 $H^{k-m}$ 进行某种估计，即可在展开深度为 $m$ 的情况下对 $d^k=-H^k\nabla f(x^k)$ 进行近似求解。当用数量矩阵来近似时，即 $\hat{H}^{k-m}=\gamma_k I$ ，其中 $\displaystyle\gamma_k=\frac{(s^{k-1})^\top y^{k-1}}{(y^{k-1})^\top y^{k-1}}$ 对应 BB 方法的第二个步长。

第一个循环：初始化 $q^{k}=-\nabla f(x^k)$ ，迭代 $\alpha_i=\rho_i(s^i)^\top q^{i+1}$ ， $q^{i}=q^{i+1}-\alpha_i y^i$ 。其中 $i$ 从 $k-1$ 减小到 $k-m$ 。不难证明 $q^{k-m}$ 递归地求出了 $-V^{k-m}\cdots V^{k-1}\nabla f(x^k)$ 并同时求出了 $\alpha_{i}=-\rho_{i}(s^{i})^\top (V^{i+1}\cdots V^{k-1})\nabla f(x^k)$ 。

经过第一个循环，我们可以将上述展开的表达式重写为

$\displaystyle \begin{array}{rl} -H^k\nabla f(x^k)= & -(V^{k-m}\cdots V^{k-1})H^{k-m}q^{k-m}\\ & -(V^{k-m+1}\cdots V^{k-1})^\top s^{k-m}\alpha_{k-m}\\ & -(V^{k-m+2}\cdots V^{k-1})^\top s^{k-m+1}\alpha_{k-m+1}\\ &-\cdots\\ &-s^{k-1}\alpha_{k-1}.\end{array}$

再引入第二个循环对这一求和式进行计算。

初始化 $r^{k-m}=\hat{H}^{k-m}q^{k-m}$ ，迭代

$\begin{array}{rl} \beta_{i}=&\hspace{-0.5em}\rho_i(y^i)^\top r^i,\\ r^{i+1}=&\hspace{-0.5em}r^i+(\alpha_i-\beta_i)s^i. \end{array}$

可以验证每一次循环 $r$ 都相当于对求和式提取公因式后的从前到后进行一步累加，最终得到 $r^k$ 即为所需的 $-H^k\nabla f(x^k)$ 。

利用上面双循环算法对下降方向进行求解的方法即为 L-BFGS 方法，这一方法只需记录 $m$ 步的信息 $\{(s^i,y^i)\}_{i=k-m}^{k}$ ，在每一步更新时，将新得到的 $(s^{k+1},y^{k+1})$ 覆盖 $(s^{k-m},y^{k-m})$ ，因此需要的空间大大减小。

该函数完成上述的 L-BFGS 算法，利用双循环算法计算下降方向，并利用线搜索确定步长。

初始化和迭代准备

函数输入： x 为迭代的初始点， fun 提供函数值和梯度， opts 为提供算法参数的结构体。

函数输出： x 为迭代得到的解， f 和 g 为该点处的函数值和梯度， Out 为记录迭代信息的结构体。

Out.f ：迭代过程的函数值信息
Out.nrmG ：迭代过程的梯度范数信息
Out.xitr ：迭代过程的优化变量 $x$ （仅在 storeitr 不为 0 时存在）
Out.nfe ：调用目标函数的次数
Out.nge ：调用梯度的次数
Out.nrmg ：迭代终止时的梯度范数
Out.iter ：迭代步数
Out.msg ：标记是否达到收敛

function [x, f, g, Out]= fminLBFGS_Loop(x, fun, opts, varargin)

从输入的结构体 opts 中读取参数或采取默认参数。

opts.xtol ：主循环针对优化变量的停机判断依据
opts.gtol ：主循环针对梯度范数的停机判断依据
opts.ftol ：主循环针对函数值的停机判断依据
opts.rho1 ：线搜索标准 $\rho_1$
opts.rho2 ：线搜索标准 $\rho_2$
opts.m ：L-BFGS 的内存对存储数
opts.maxit ：主循环的最大迭代次数
opts.storeitr ：标记是否记录每一步迭代的 $x$
opts.record ：标记是否需要迭代信息的输出
opts.itPrint ：每隔几步输出一次迭代信息

if ~isfield(opts, 'xtol');      opts.xtol = 1e-6; end
if ~isfield(opts, 'gtol');      opts.gtol = 1e-6; end
if ~isfield(opts, 'ftol');      opts.ftol = 1e-16; end
if ~isfield(opts, 'rho1');      opts.rho1  = 1e-4; end
if ~isfield(opts, 'rho2');      opts.rho2  = 0.9; end
if ~isfield(opts, 'm');         opts.m  = 5; end
if ~isfield(opts, 'maxit');     opts.maxit  = 1000; end
if ~isfield(opts, 'storeitr');  opts.storeitr = 0; end
if ~isfield(opts, 'record');    opts.record = 0; end
if ~isfield(opts,'itPrint');    opts.itPrint = 1;   end

参数复制。

xtol = opts.xtol;
ftol = opts.ftol;
gtol = opts.gtol;
maxit = opts.maxit;
storeitr = opts.storeitr;
parsls.ftol = opts.rho1;
parsls.gtol = opts.rho2;
m = opts.m;
record = opts.record;
itPrint = opts.itPrint;

初始化和迭代准备，计算初始点处的信息。初始化迭代信息。

[f,  g] = feval(fun, x , varargin{:});
nrmx = norm(x);
Out.f = f;  Out.nfe = 1; Out.nrmG = [];

在 storeitr 不为 0 时，记录每一步迭代的 $x$ 。

if storeitr
    Out.xitr = x;
end

SK , YK 用于存储 L-BFGS 算法中最近的 $m$ 步的 $s$ （ $x$ 的变化量）和 $y$ （梯度 $g$ 的变化量）。

n = length(x);
SK = zeros(n,m);
YK = zeros(n,m);
istore = 0; pos = 0;  status = 0;  perm = [];

为打印每一步的迭代信息设定格式。

if record == 1
    if ispc; str1 = '  %10s'; str2 = '  %6s';
    else     str1 = '  %10s'; str2 = '  %6s'; end
    stra = ['%5s',str2,str2,str1, str2, str2,'\n'];
    str_head = sprintf(stra, ...
        'iter', 'stp', 'obj', 'diffx', 'nrmG', 'task');
    str_num = ['%4d  %+2.1e  %+2.1e  %+2.1e  %+2.1e  %2d\n'];
end

迭代主循环

迭代最大步数 maxit 。当未达到收敛条件时，记录为超过最大迭代步数退出。

Out.msg = 'MaxIter';
for iter = 1:maxit

记录上一步迭代的结果。

    xp = x;   nrmxp = nrmx;
    fp = f;   gp = g;

L-BFGS 双循环方法寻找下降方向。在第一次迭代时采用负梯度方向，之后便使用 L-BFGS 方法来估计 $d=-Hg$ 。

    if istore == 0
        d = -g;
    else
        d = LBFGS_Hg_Loop(-g);
    end

沿 L-BFGS 方法得到的下降方向做线搜索。调用函数 ls_csrch 进行线搜索，其参考了 MINPACK-2 中的线搜索函数。

首先初始化线搜索标记 workls.task 为 1， deriv 为目标函数沿当前下降方向的方向导数。通过线搜索寻找合适的步长 $\alpha_k$ ，使得以下条件满足：

$\begin{array}{rl} f(x^k+\alpha_k d^k)&\hspace{-0.5em}\le f(x^k)+\rho_1\cdot\alpha_kg(x^k), \\ |g(x^k+\alpha_kd^k)|&\hspace{-0.5em}\le \rho_2\cdot |g(x^k)|. \end{array}$

ls_csrch 每次调用只执行线搜索的一步，并用 workls.task 指示下一步应当执行的操作。此处 workls.task==2 意味着需要重新计算当前点函数值和梯度等。具体的步长寻找过程比较复杂，可以参考相应文件。

直到满足线搜索条件时，退出线搜索循环，得到更新之后的 $x$ 。

    workls.task =1;
    deriv = d'*g;
    normd = norm(d);

    stp = 1;
    while 1
        [stp, f, deriv, parsls, workls] = ....
            ls_csrch(stp, f, deriv , parsls , workls);

        if (workls.task == 2)
            x = xp + stp*d;
            [f,  g] = feval(fun, x, varargin{:});
            Out.nfe = Out.nfe + 1;
            deriv = g'*d;
        else
            break
        end
    end

对于线搜索得到的步长 $\alpha_k$ ，令 $s^k=x^{k+1}-x^k=\alpha_kd^k$ ，则 $\|s^k\|_2=\alpha_k\|d^k\|_2$ 。计算 $\|s^k\|_2/\max(1,\|x^k\|_2)$ 并将其作为判断收敛的标准。

    nrms = stp*normd;
    diffX = nrms/max(nrmxp,1);

更新 $\|x^{k+1}\|_2$ , $\|g^{k+1}\|_2$ ，记录一步迭代信息。

    nrmG =  norm(g);
    Out.nrmg =  nrmG;
    Out.f = [Out.f; f];
    Out.nrmG = [Out.nrmG; nrmG];
    if storeitr
        Out.xitr = [Out.xitr, x];
    end
    nrmx = norm(x);

停机准则， diffX 表示相邻迭代步 $x$ 的相对变化， nrmG 表示当前 $x$ 处的梯度范数， $\displaystyle\frac{|f(x^{k+1})-f(x^k)|}{1+|f(x^k)|}$ 用以表示函数值的相对变化。当前两者均小于阈值，或者第三者小于阈值时，认为达到停机标准，退出当前循环。

    cstop = ((diffX < xtol) && (nrmG < gtol) )|| (abs(fp-f)/(abs(fp)+1)) < ftol;

当需要详细输出时，在（1）开始迭代时（2）达到收敛时（3）达到最大迭代次数或退出迭代时（4）每若干步，打印详细结果。

    if (record == 1) && (cstop || iter == 1 || iter==maxit || mod(iter,itPrint)==0)
        if iter == 1 || mod(iter,20*itPrint) == 0 && iter~=maxit && ~cstop
            fprintf('\n%s', str_head);
        end
        fprintf(str_num, ...
            iter, stp, f, diffX, nrmG, workls.task);
    end

当达到收敛条件时，停止迭代，记为达到收敛。

    if cstop
        Out.msg = 'Converge';
        break;
    end

计算 $s^k=x^{k+1}-x^{k}$ , $y^k=g^{k+1}-g^{k}$ 。当得到的 $\|y^k\|$ 不小于阈值时，保存当前的 $s^k, y^k$ ，否则略去。利用 pos 记录当前存储位置，然后覆盖该位置上原来的信息。

    ygk = g-gp;		s = x-xp;
    if ygk'*ygk>1e-20
        istore = istore + 1;
        pos = mod(istore, m); if pos == 0; pos = m; end
        YK(:,pos) = ygk;  SK(:,pos) = s;   rho(pos) = 1/(ygk'*s);

用于提供给 L-BFGS 双循环递归算法，以指明双循环的循环次数。当已有的记录超过 $m$ 时，则循环 $m$ 次。否则，循环次数等于当前的记录个数。 perm 按照顺序记录存储位置。

        if istore <= m; status = istore; perm = [perm, pos];
        else status = m; perm = [perm(2:m), perm(1)]; end
    end
end

当从上述循环中退出时，记录输出。

Out.iter = iter;
Out.nge = Out.nfe;

L-BFGS 双循环递归算法

利用双循环递归算法，返回下一步的搜索方向即 $-Hg$ 。初始化 $q$ 为初始方向，在 L-BFGS 主算法中，这一方向为负梯度方向。

    function y = LBFGS_Hg_Loop(dv)
        q = dv;   alpha = zeros(status,1);

第一个循环， status 步迭代。（ status 的计算见上，当迭代步足够大时为 $m$ ） perm 按照顺序记录了存储位置。从中提取出位置 $k$ 的格式为： $\alpha_i=\rho_i(s^i)^\top q^{i+1}$ , $q^{i}=q^{i+1}-\alpha_i y^i$ 。其中 $i$ 从 $k-1$ 减小到 $k-m$ 。

        for di = status:-1:1
            k = perm(di);
            alpha(di) = (q'*SK(:,k)) * rho(k);
            q = q - alpha(di)*YK(:,k);
        end

$r^{k-m}=\hat{H}^{k-m}q^{k-m}.$

        y = q/(rho(pos)* (ygk'*ygk));

第二个循环，迭代格式 $\beta_{i}=\rho_i(y^i)^\top r^i$ , $r^{i+1}=r^i+(\alpha_i-\beta_i)s^i$ 。代码中的 y 对应于迭代格式中的 $r$ ，当两次循环结束时，以返回的 y 的值作为下降方向。

        for di = 1:status
            k = perm(di);
            beta = rho(k)* (y'* YK(:,k));
            y = y + SK(:,k)*(alpha(di)-beta);
        end
    end

end

参考页面

关于 L-BFGS 算法的应用，参考页面实例：L-BFGS 解逻辑回归问题以及实例：L-BFGS 解基追踪问题。

我们在其中利用了翻译为 MATLAB 代码的 MINPACK-2 的线搜索函数 ls_csrch ，代码详见线搜索函数和被其调用的线搜索辅助函数，或参考 Fortran 版本的官方代码 MINPACK-2。

此页面的源代码请见： fminLBFGS_Loop.m。

版权声明