LASSO 问题的分块坐标下降法

利用分块坐标下降法（BCD）优化如下的 LASSO 问题

$\displaystyle\min_x \frac{1}{2}\|Ax-b\|_2^2 + \mu \|x\|_1.$

该算法被外层连续化策略函数调用，完成对某一固定正则化系数 $\mu_k$ 的内层迭代优化。

分块坐标下降法将 $x$ 的第 $i$ 个分量 $x_i$ 作为第 $i$ 块变量。记 $a_i$ 为 $A$ 的第 $i$ 列， $\bar{x}_i$ 和 $\bar{A}_i$ 为 $x$ 去掉第 $i$ 个分量和 $A$ 去掉第 $i$ 列后的量，则在第 $i$ 块的更新中，原问题等价于

$\displaystyle\min_{x_i}\mu|x_i|+\mu\|\bar{x}_i\|_1 +\frac{1}{2}\|a_ix_i-(b-\bar{A}_i\bar{x}_i)\|_2^2.$

记 $c_i=b-\bar{A}_i\bar{x}_i$ ，则由上述问题可以得到迭代格式

$\displaystyle x_i^k= \left\{\begin{array}{ll} \displaystyle\frac{a_i^\top c_i-\mu}{\|a_i\|^2_2}, & a_i^\top c_i>\mu;\\ \displaystyle\frac{a_i^\top c_i+\mu}{\|a_i\|^2_2}, & a_i^\top c_i<-\mu;\\ 0, & \mathrm{otherwise}. \end{array}\right.$

初始化和迭代准备

函数在 LASSO 连续化策略下，完成内层迭代的优化。

输入信息： $A$ , $b$ ，当前内层迭代的正则化系数 $\mu$ ，迭代初始值 $x^0$ ，原问题对应的正则化系数 $\mu_0$ ，以及提供各参数的结构体 opts 。

输出信息：迭代得到的解 $x$ 和结构体 out 。

out.fvec ：每一步迭代的原始目标函数值（对应于原问题的 $\mu_0$ ）
out.fval ：迭代终止时的原始目标函数值（对应于原问题的 $\mu_0$ ）
out.nrmG ：迭代终止时的梯度范数
out.tt ：运行时间
out.itr ：迭代次数
out.flag ：标记是否达到收敛

function [x, out] = LASSO_bcd_inn(x0, A, b, mu, mu0, opts)

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

opts.maxit ：最大迭代次数
opts.ftol ：针对函数值的停机准则，当相邻两次迭代函数值之差小于该值时认为该条件满足
opts.gtol ：针对梯度的停机准则，当当前步梯度范数小于该值时认为该条件满足
opts.verbose ：不为 0 时输出每步迭代信息，否则不输出

if ~isfield(opts, 'maxit'); opts.maxit = 200; 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

迭代准备，计算初始时刻 $x$ 对应的各变量值。

out = struct();
x = x0;
tt = tic;
[m,n] = size(A);
Ax = A*x;
r = Ax - b;
f = .5*norm(r,2)^2 + mu0*norm(x,1);

记录初始时刻的目标函数值（对应原正则化系数 mu_0）。

out.fvec = f;

迭代主循环

内层迭代，以 maxit 为最大迭代次数。

for k = 1:opts.maxit

记录上一步的迭代结果。

    fp = f;

BCD 的迭代过程，每一步迭代分别针对 $x$ 的第 $i$ 个坐标分量进行优化，计算 $c_i=b-\bar{A}_i\bar{x}_i$ 。

    for i = 1:n
        a = A(:,i);
        Axnew = Ax - a*x(i);
        c = b - Axnew;

BCD 的每一个坐标块的下降。

$\displaystyle x_i^k= \left\{\begin{array}{ll} \displaystyle\frac{a_i^\top c_i-\mu}{\|a_i\|^2_2}, & a_i^\top c_i>\mu;\\ \displaystyle\frac{a_i^\top c_i+\mu}{\|a_i\|^2_5}, & a_i^\top c_i<-\mu;\\ 0, & \mathrm{otherwise}. \end{array}\right.$

        ac = a'*c;
        if ac > mu
            x(i) = (ac - mu)/norm(a,2)^2;
        elseif ac < -mu
            x(i) = (ac + mu)/norm(a,2)^2;
        else
            x(i) = 0;
        end

$Ax=\bar{A}_i\bar{x}_i+a_ix_i.$

        Ax = Axnew + a*x(i);
    end

从上面的循环中退出， $x$ 的全部 $n$ 个坐标分量依次优化完成，更新并记录相关变量的值。

    Axb = Ax - b;
    f = .5*norm(Axb,2)^2 + mu0*norm(x,1);

利用一步近似点梯度法作为梯度的近似。

    nrmG = norm(x - prox(x - A'*Axb, mu));

记录一步 BCD 迭代的函数值。

    out.fvec = [out.fvec, f];

当 verbose 不为 0 时，输出迭代信息。

    if opts.verbose
        fprintf('itr: %d\t fval: %e \t nrmG: %.1e \n', k, f, nrmG);
    end

停机准则，当函数值变化小于阈值或者梯度范数小于阈值时，认为达到收敛，退出迭代循环。当退出循环时，向外层迭代（连续化策略）报告内层迭代的退出方式，当达到最大迭代次数退出时， out.flag 记为 1，否则则为达到收敛，记为 0。这个指标用于判断是否进行正则化系数的衰减。

    if abs(f-fp) < opts.ftol || nrmG < opts.gtol
        break;
    end
end

if k == opts.maxit
    out.flag = 1;
else
    out.flag = 0;
end

记录内层迭代的输出。

out.fvec = out.fvec(1:k);
out.fval = f;
out.itr = k;
out.tt = toc(tt);

end

辅助函数

函数 $h(x)=\mu\|x\|_1$ 对应的邻近算子 $\mathrm{sign}(x)\max\{|x|-\mu,0\}$ 。

function y = prox(x, mu)
y = max(abs(x) - mu, 0);
y = sign(x) .* y;
end

参考页面

该函数由连续化策略调用，关于连续化策略参见 LASSO问题连续化策略。

我们在页面实例：利用分块坐标下降法求解 LASSO 问题中构建一个 LASSO 问题并展示该算法的一个应用。

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

版权声明