LASSO 问题的近似点算法
近似点算法(PPA)对于 LASSO 问题,考虑其等价形式
近似点算法的一个迭代步为
其中 。对于子问题考虑其对偶问题,引入拉格朗日乘子
,并令
则其迭代格式满足
以 为逼近解,
和
的子问题是半光滑的, 因此利用半光滑牛顿法加速的梯度法进行求解。 再结合最优性条件,有
目录
初始化和迭代准备
PPA不使用连续化策略,直接对原始的 LASSO 问题进行求解。
输入信息: ,
,
,迭代初始值
以及提供各参数的结构体 opts 。
输出信息: 迭代得到的解 和结构体 out 。
- out.fvec :每一步迭代的 LASSO 目标函数值
- out.flag :标记是否达到收敛
- out.fval :迭代终止时的 LASSO 目标函数值
- out.tt :运行时间
- out.itr :外层迭代次数
function [x, out] = LASSO_ppa(x0, A, b, mu, opts)
从输入的结构体 opts 中读取参数或采取默认参数。
- opts.maxit :最大迭代次数
- opts.ftol :针对函数值的停机准则,当相邻两次迭代函数值之差小于该值时认为该条件满足
- opts.gtol :针对梯度的停机准则,当当前步梯度范数小于该值时认为该条件满足
- opts.verbose : >1 时输出每一步信息, 0-1 时输出外层迭代信息, =0 时不输出
- opts.t0 :初始步长
if ~isfield(opts, 'maxit'); opts.maxit = 500; end if ~isfield(opts, 'ftol'); opts.ftol = 1e-8; end if ~isfield(opts, 'gtol'); opts.gtol = 1e-6; end if ~isfield(opts, 'verbose'); opts.verbose = 0; end if ~isfield(opts, 't0'); opts.t0 = 1e3; end out = struct();
optsz 为 z 的子问题提供参数。
optsz = struct(); optsz.verbose = (opts.verbose > 1);
迭代初始化。
以 x^0 为迭代初始点,在初始 x^0 处计算相关变量的值。
x = x0; fp = inf; tt = tic; r = A*x - b; f = .5*norm(r,2)^2 + mu*norm(x,1); out.fvec = f;
t 为步长,初始值设定为 opts.t0。
t = opts.t0;
变量 y=Ax-b 和拉格朗日乘子 z。
z = zeros(size(b)); y = A*x - b; toldiff = 1;
迭代主循环
以 maxit 为最大迭代步数。每次迭代开始,记录上一次迭代的目标函数值。
for k = 1:opts.maxit
fp = f;
toldiff 为 ,这里为内层的
的子问题迭代确定停机准则:
,并令其中的
,
。
optsz.gtol = 8/k^2*min(1, toldiff);
内层循环的梯度法解拉格朗日乘子 的子问题。然后对
,
更新。
[z,~] = zgbb(z,A,b,mu,t, x,y, optsz); xp = x; x = prox(x - t*A'*z, mu*t); yp = y; y = 1/(t+1)*(y + t*z);
一步 PPA 的外层迭代结束,计算目标函数值,并以一步近似点梯度法的 估计梯度。
Axb = A*x - b; f = .5*norm(Axb,2)^2 + mu*norm(x,1); nrmG = norm(x - prox(x - A'*Axb, mu)); toldiff = norm(x-xp,2) + norm(y- yp,2);
记录一步迭代的目标函数值。
out.fvec = [out.fvec; f];
当 opts.verbose 不为 0 时输出详细的迭代信息。
if opts.verbose fprintf('itr: %d\t t: %e\t fval: %e \t nrmG: %.1e \t feasi: %.1e \n', ... k, t, f, nrmG, norm(A*x -b - y,2)); end
外层迭代的停机准则,当函数值的变化或梯度(一步近似点梯度法估计)范数小于阈值时, 认为收敛,退出循环。
if abs(f-fp) < opts.ftol || nrmG < opts.gtol break; end end
当收敛时,记录输出信息。以 out.flag 记录迭代结束原因, out.flag=1 表示满足梯度条件, out.flag=0 表示满足函数值条件。
if abs(f - fp) < opts.ftol out.flag = 0; else out.flag = 1; end
记录输出信息。
out.fval = f; out.itr = k; out.tt = toc(tt);
end
拉格朗日乘子
子问题
从输入的结构体中读取参数或采取默认参数。
- optsz.maxit :最大迭代次数
- optsz.gtol :子问题的停机准则,当梯度范数小于该值时认为满足; 注意到该值随着外层迭代而变化,具体的计算方法见上文
- optsz.verbose :不为 0 时输出每步迭代信息,否则不输出
- optsz.tau0 :初始步长
function [z,out] = zgbb(z0,A,b,mu,t,x,y, optsz) if ~isfield(optsz, 'maxit'); optsz.maxit = 500; end if ~isfield(optsz, 'gtol'); optsz.gtol = 1e-6; end if ~isfield(optsz, 'verbose'); optsz.verbose = 0; end if ~isfield(optsz, 'tau0'); optsz.tau0 = 1e-2; end
初始化变量。
out = struct(); [m,n] = size(A); Im = eye(m); H = eye(m); z = z0; r = x - t*(A'*z); w = prox(r,mu*t);
初始步长设定为 optsz.tau0。
tau = optsz.tau0;
将子问题中的 的部分解析地代入,得到其迭代的目标函数(去掉常数项)
其中 。
令 ,
。
f = - mu*(norm(w,1) + 1/(2*mu*t)*norm(w-r,2)^2) + 1/(2*t)*norm(r,2)^2 ...
+ t/(2*(t+1))*(norm(z,2)^2) + 1/(t+1)*(y'*z) + b'*z;
初始化线搜索参数。
Q = 1; Cval = f; gamma = 0.85; rhols = 1e-6;
目标函数为连续可微函数,且其梯度为 。
g = - A*w + t/(t+1)*z + 1/(t+1)*y + b; nrmg = norm(g,2);
对 的梯度法迭代主循环,最大迭代步 maxit。
maxit = optsz.maxit;
for k = 1:maxit
zp = z;
gp = g;
进行线搜索。 nls 记录线搜索次数, 表示每次不满足线搜索时步长减小的比例。 每次循环以当前步长进行一步试探,下降方向
,求出在新的
处的变量和目标函数值。 当
时,则为一般的带线搜索的 BB 步长梯度法,从 10 步后,
为广义海瑟矩阵的逆,步长设定为
,对应半光滑牛顿法。
nls = 1; deriv = rhols*nrmg^2; eta = 0.2; while 1 z = zp - tau*(H*g); r = x - t*(A'*z); w = prox(r, mu*t); f = - mu*(norm(w,1) + 1/(2*mu*t)*norm(w-r,2)^2) + 1/(2*t)*norm(r,2)^2 ... + t/(2*(t+1))*(norm(z,2)^2) + 1/(t+1)*(y'*z) + b'*z;
满足线搜索准则 (Zhang & Hager) 或进行超过 10 次步长衰减后退出线搜索,否则以
的比例对步长进行衰减。
if f < Cval - tau*deriv || nls == 10 break; end tau = eta*tau; nls = nls + 1; end
更新梯度。
g = - A*w + t/(t+1)*z + 1/(t+1)*y + b; nrmg = norm(g,2);
当 optsz.verbose 不为 0 时输出详细的迭代信息。
if optsz.verbose fprintf('iter: %d \t f: %.4e \t nrmG: %.1e \t nls: %d \n', k, f, nrmg, nls); end
停机准则(见上文),当满足停机准则时退出循环。
gtol = sqrt(1/(t+1))*optsz.gtol; if nrmg < gtol break; end
BB 步长的计算,令 ,
,则两种 BB 步长
,
。 这里我们交替使用两种 BB 步长并限定在
中。
dz = z - zp; dg = g - gp; dzg = abs(dz'*dg); if dzg > 0 if mod(k,2) == 0 tau = norm(dz,2)^2 /dzg; else tau = dzg/(norm(dg,2)^2); end end tau = min(max(1e-12,tau),1e12);
线搜索参数的更新。
Qp = Q; Q = gamma*Qp + 1; Cval = (gamma*Qp*Cval + f)/Q;
从 10 步之后,开始使用半光滑牛顿法。 由优化的目标函数的形式
得到其广义海瑟矩阵
其中 为满足
的下标对应的
的列。再由 Sherman-Morrison-Woodbury 公式可得
为广义海瑟矩阵的逆。在利用半光滑牛顿法的情况下,采用 的固定步长和一般的线搜索准则 (即 Cval 直接取当前点函数值)。
if k > 10 ind = find(abs(r) > mu*t); Aind = A(:,ind); tmp = 1/(t+1)*eye(length(ind)) + Aind'*Aind; H = (t+1)/t*(Im - Aind*(tmp\Aind')); tau = 1; Cval = f; end end
当从迭代中退出时,返回这一内重循环的迭代步。
out.iter = k;
end
辅助函数
函数 对应的邻近算子
。
function y = prox(x, mu) y = max(abs(x) - mu, 0); y = sign(x) .* y; end
参考页面
在页面 实例:近似点法解 LASSO 问题 中我们构建一个 LASSO 问题并展示该算法在其中的应用。
此页面的源代码请见: LASSO_ppa.m。
版权声明
此页面为《最优化:建模、算法与理论》、《最优化计算方法》配套代码。 代码作者:文再文、刘浩洋、户将,代码整理与页面制作:杨昊桐。
著作权所有 (C) 2020 文再文、刘浩洋、户将