L-BFGS算法基础_matlab
L-BFGS 优化基础
针对无约束优化问题： min x f ( x ) , \min_x f(x), xminf(x),
L-BFGS 在拟牛顿法 BFGS 迭代格式的基础上进行修改，用以解决大规模问题的存储和计算困难。对于拟牛顿法中的迭代方向 d k = − H k ∇ f ( x k ) d^k=-H^k\nabla f(x^k) dk=−Hk∇f(xk)。
考虑利用递归展开的方式进行求解
首先，对于 BFGS 迭代格式， H k + 1 = ( V k ) ⊤ H k V k + ρ k s k ( s k ) ⊤ H^{k+1}=(V^k)^\top H^kV^k+\rho_ks^k(s^k)^\top Hk+1=(Vk)⊤HkVk+ρksk(sk)⊤，其中 ρ k = 1 ( y k ) ⊤ s k , \displaystyle\rho_k = \frac{1}{(y^k)^\top s^k}, ρk=(yk)⊤sk1, V k = I − ρ k y k ( s k ) ⊤ V^k=I-\rho_k y^k(s^k)^\top Vk=I−ρkyk(sk)⊤。将其递归地展开得到：
− H k ∇ f ( x k ) = − ( V k − m ⋯ V k − 1 ) ⊤ H k − m ( V k − m ⋯ V k − 1 ) ∇ f ( x k ) − ρ k − m ( V k − m + 1 ⋯ V k − 1 ) ⊤ s k − m ( s k − m ) ⊤ ( V k − m + 1 ⋯ V k − 1 ) ∇ f ( x k ) − ρ k − m + 1 ( V k − m + 2 ⋯ V k − 1 ) ⊤ s k − m + 1 ( s k − m + 1 ) ⊤ ( V k − m + 2 ⋯ V k − 1 ) ∇ f ( x k ) − ⋯ − ρ k − 1 s k − 1 ( s k − 1 ⊤ ) ∇ f ( x k ) . \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} −Hk∇f(xk)=−(Vk−m⋯Vk−1)⊤Hk−m(Vk−m⋯Vk−1)∇f(xk)−ρk−m(Vk−m+1⋯Vk−1)⊤sk−m(sk−m)⊤(Vk−m+1⋯Vk−1)∇f(xk)−ρk−m+1(Vk−m+2⋯Vk−1)⊤sk−m+1(sk−m+1)⊤(Vk−m+2⋯Vk−1)∇f(xk)−⋯−ρk−1sk−1(sk−1⊤)∇f(xk).
我们只需对其中的 H k − m H^{k-m} Hk−m 进行某种估计，即可在展开深度为 m m m 的情况下对 d k = − H k ∇ f ( x k ) d^k=-H^k\nabla f(x^k) dk=−Hk∇f(xk) 进行近似求解。当用数量矩阵来近似时，即 H ^ k − m = γ k I \hat{H}^{k-m}=\gamma_k I H^k−m=γkI，其中 γ k = ( s k − 1 ) ⊤ y k − 1 ( y k − 1 ) ⊤ y k − 1 \displaystyle\gamma_k=\frac{(s^{k-1})^\top y^{k-1}}{(y^{k-1})^\top y^{k-1}} γk=(yk−1)⊤yk−1(sk−1)⊤yk−1 对应 BB 方法的第二个步长。
第一个循环：
初始化
q
k
=
−
∇
f
(
x
k
)
q^{k}=-\nabla f(x^k)
qk=−∇f(xk)，迭代
α
i
=
ρ
i
(
s
i
)
⊤
q
i
+
1
\alpha_i=\rho_i(s^i)^\top q^{i+1}
αi=ρi(si)⊤qi+1，
q
i
=
q
i
+
1
−
α
i
y
i
q^{i}=q^{i+1}-\alpha_i y^i
qi=qi+1−αiyi。其中
i
i
i 从
k
−
1
k-1
k−1 减小到
k
−
m
k-m
k−m。
不难证明
q
k
−
m
=
−
V
k
−
m
⋯
V
k
−
1
∇
f
(
x
k
)
q^{k-m}=-V^{k-m}\cdots V^{k-1}\nabla f(x^k)
qk−m=−Vk−m⋯Vk−1∇f(xk) 并同时求出了
α
i
=
−
ρ
i
(
s
i
)
⊤
(
V
i
+
1
⋯
V
k
−
1
)
∇
f
(
x
k
)
\alpha_{i}=-\rho_{i}(s^{i})^\top (V^{i+1}\cdots V^{k-1})\nabla f(x^k)
αi=−ρi(si)⊤(Vi+1⋯Vk−1)∇f(xk)。
经过第一个循环，我们可以将上述展开的表达式重写为：
− H k ∇ f ( x k ) = − ( V k − m ⋯ V k − 1 ) H k − m q k − m − ( V k − m + 1 ⋯ V k − 1 ) ⊤ s k − m α k − m − ( V k − m + 2 ⋯ V k − 1 ) ⊤ s k − m + 1 α k − m + 1 − ⋯ − s k − 1 α k − 1 . \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} −Hk∇f(xk)=−(Vk−m⋯Vk−1)Hk−mqk−m−(Vk−m+1⋯Vk−1)⊤sk−mαk−m−(Vk−m+2⋯Vk−1)⊤sk−m+1αk−m+1−⋯−sk−1αk−1.
第二个循环：
对这一求和式进行计算。
初始化 r k − m = H ^ k − m q k − m r^{k-m}=\hat{H}^{k-m}q^{k-m} rk−m=H^k−mqk−m，迭代： β i = ρ i ( y i ) ⊤ r i , r i + 1 = r i + ( α i − β i ) s i . \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} βi=ri+1=ρi(yi)⊤ri,ri+(αi−βi)si.
可以验证每一次循环
r
r
r 都相当于对求和式提取公因式后的从前到后进行一步累加，最终得到
r
k
r^k
rk 即为所需的
−
H
k
∇
f
(
x
k
)
-H^k\nabla f(x^k)
−Hk∇f(xk)。
利用上面双循环算法对下降方向进行求解的方法即为 L-BFGS 方法，这一方法只需记录 m m m 步的信息 { ( s i , y i ) } i = k − m k \{(s^i,y^i)\}_{i=k-m}^{k} {(si,yi)}i=k−mk，在每一步更新时，将新得到的 ( s k + 1 , y k + 1 ) (s^{k+1},y^{k+1}) (sk+1,yk+1) 覆盖 ( s k − m , y k − m ) (s^{k-m},y^{k-m}) (sk−m,yk−m)，因此需要的空间大大减小。
L-BFGS 算法_matlab
该函数完成上述的 L-BFGS 算法，利用双循环算法计算下降方向，并利用线搜索确定步长。
初始化和迭代准备
函数输入： x 为迭代的初始点， fun 提供函数值和梯度， opts 为提供算法参数的结构体。
函数输出： x 为迭代得到的解， f 和 g 为该点处的函数值和梯度， Out
为记录迭代信息的结构体。
- Out.f ：迭代过程的函数值信息
- Out.nrmG ：迭代过程的梯度范数信息
- Out.xitr ：迭代过程的优化变量 x x 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 ：线搜索标准 ρ 1 \rho_1 ρ1
- opts.rho2 ：线搜索标准 ρ 2 \rho_2 ρ2
- opts.m ：L-BFGS 的内存对存储数
- opts.maxit ：主循环的最大迭代次数
- opts.storeitr ：标记是否记录每一步迭代的 x x 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); # x的1范数
Out.f = f;
Out.nfe = 1;
Out.nrmG = [];
在 storeitr 不为 0 时，记录每一步迭代的 x x x。
if storeitr
Out.xitr = x;
end
SK , YK 用于存储 L-BFGS 算法中最近的 m m m 步的 s s s（ x x x 的变化量）和 y y y（梯度 g g 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
=
−
H
g
d=-Hg
d=−Hg。
if istore == 0
d = -g;
else
d = LBFGS_Hg_Loop(-g);
end
沿 L-BFGS 方法得到的下降方向做线搜索。调用函数 ls_csrch 进行线搜索，其参考了MINPACK-2 中的线搜索函数。
首先初始化线搜索标记 workls.task 为 1， deriv 为目标函数沿当前下降方向的方向导数。通过线搜索寻找合适的步长 α k \alpha_k αk，使得以下条件满足： f ( x k + α k d k ) ≤ f ( x k ) + ρ 1 ⋅ α k g ( x k ) , ∣ g ( x k + α k d k ) ∣ ≤ ρ 2 ⋅ ∣ g ( x 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} f(xk+αkdk)∣g(xk+αkdk)∣≤f(xk)+ρ1⋅αkg(xk),≤ρ2⋅∣g(xk)∣.
ls_csrch 每次调用只执行线搜索的一步，并用 workls.task 指示下一步应当执行的操作。
此处 workls.task==2 意味着需要重新计算当前点函数值和梯度等。具体的步长寻找过程比较复杂，可以参考相应文件。
直到满足线搜索条件时，退出线搜索循环，得到更新之后的 x x 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
对于线搜索得到的步长 α k \alpha_k αk，令 s k = x k + 1 − x k = α k d k s^k=x^{k+1}-x^k=\alpha_kd^k sk=xk+1−xk=αkdk，则 ∥ s k ∥ 2 = α k ∥ d k ∥ 2 \|s^k\|_2=\alpha_k\|d^k\|_2 ∥sk∥2=αk∥dk∥2。计算 ∥ s k ∥ 2 / max ( 1 , ∥ x k ∥ 2 ) \|s^k\|_2/\max(1,\|x^k\|_2) ∥sk∥2/max(1,∥xk∥2)，并将其作为判断收敛的标准。
nrms = stp*normd;
diffX = nrms/max(nrmxp,1);
更新 ∥ x k + 1 ∥ 2 \|x^{k+1}\|_2 ∥xk+1∥2, ∥ g k + 1 ∥ 2 \|g^{k+1}\|_2 ∥gk+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 x x 的相对变化， nrmG 表示当前 x x x 处的梯度范数， ∣ f ( x k + 1 ) − f ( x k ) ∣ 1 + ∣ f ( x k ) ∣ \displaystyle\frac{|f(x^{k+1})-f(x^k)|}{1+|f(x^k)|} 1+∣f(xk)∣∣f(xk+1)−f(xk)∣用以表示函数值的相对变化。当前两者均小于阈值，或者第三者小于阈值时，认为达到停机标准，退出当前循环。
cstop = ((diffX < xtol) && (nrmG < gtol) )|| (abs(fp-f)/(abs(fp)+1)) < ftol;
当需要详细输出时，在开始迭代时、达到收敛时、达到最大迭代次数或退出迭代时、每若干步，打印详细结果。
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
s^k=x^{k+1}-x^{k}
sk=xk+1−xk,
y
k
=
g
k
+
1
−
g
k
y^k=g^{k+1}-g^{k}
yk=gk+1−gk。当得到的
∥
y
k
∥
\|y^k\|
∥yk∥ 不小于阈值时，保存当前的
s
k
,
y
k
s^k, y^k
sk,yk，否则略去。利用 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 m 时，则循环 m m m 次。否则，循环次数等于当前的记录个数。 perm 按照顺序记录存储位置。
perm = [perm(2:m), perm(1)];
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 双循环递归算法
利用双循环递归算法，返回下一步的搜索方向即 − H g -Hg −Hg。初始化 q q q 为初始方向，在 L-BFGS 主算法中，这一方向为负梯度方向。
function y = LBFGS_Hg_Loop(dv)
q = dv;
alpha = zeros(status,1);
第一个循环： status 步迭代。（ status 的计算见上，当迭代步足够大时为 m m m ）perm 按照顺序记录了存储位置。从中提取出位置 k k k 的格式为： α i = ρ i ( s i ) ⊤ q i + 1 \alpha_i=\rho_i(s^i)^\top q^{i+1} αi=ρi(si)⊤qi+1 , q i = q i + 1 − α i y i q^{i}=q^{i+1}-\alpha_i y^i qi=qi+1−αiyi。其中 i i i 从 k − 1 k-1 k−1 减小到 k − m k-m 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 = H ^ k − m q k − m r^{k-m}=\hat{H}^{k-m}q^{k-m} rk−m=H^k−mqk−m
y = q/(rho(pos)* (ygk’*ygk));
y = q/(rho(pos)* (ygk'*ygk));
第二个循环：迭代格式 β i = ρ i ( y i ) ⊤ r i \beta_{i}=\rho_i(y^i)^\top r^i βi=ρi(yi)⊤ri, r i + 1 = r i + ( α i − β i ) s i r^{i+1}=r^i+(\alpha_i-\beta_i)s^i ri+1=ri+(αi−βi)si。代码中的 y 对应于迭代格式中的 r r 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
440
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
