利用交替方向乘子法 (ADMM) 求解 LASSO 对偶问题

交替方向乘子法，即 the Alternating Direction Method of Multipliers （ADMM），利用 ADMM 求解 LASSO 问题的对偶问题。

针对 LASSO 问题

$\min_{x,z} \frac{1}{2} \|Ax-b\|_2^2 + \mu \|x\|_1,$

考虑其对偶问题的 ADMM 等价形式

$\begin{array}{rl} \displaystyle\min_y &\hspace{-0.5em} b^\top y +\frac{1}{2}\|y\|_2^2+I_{\|z\|_\infty\le\mu}(z),\\ \displaystyle\mathrm{s.t.} &\hspace{-0.5em} A^\top y + z = 0. \end{array}$

引入 $x$ 作为拉格朗日乘子，得到增广拉格朗日函数 $L_\rho(y,z,x)=b^\top y+\frac{1}{2}\|y\|^2+I_{\|z\|_\infty\le\mu}(z) -x^\top(A^\top y+z)+\frac{\rho}{2}\|A^\top y+z\|^2$ 。在 ADMM 的每一步迭代中，交替更新 $y$ , $z$ ，在更新 $y$ ( $z$ ) 的时候 $z$ ( $y$ ) 固定（看成常量）。

初始化和迭代准备

函数通过优化上面给出的增广拉格朗日函数，以得到 LASSO 问题的解。

输入信息： $A$ , $b$ , $\mu$ ，迭代初始值 $x^0$ 以及提供各参数的结构体 opts .

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

out.fvec ：每一步迭代的 LASSO 问题目标函数值
out.fval ：迭代终止时的 LASSO 问题目标函数值
out.tt ：运行时间
out.itr ：迭代次数
out.y ：迭代终止时的对偶变量 $y$ 的值
out.nrmC ：迭代终止时的约束违反度，在一定程度上反映收敛性

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

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

opts.maxit ：最大迭代次数
opts.ftol ：针对函数值的停机准则，当相邻两次迭代函数值之差小于该值时认为该条件满足
opts.gtol ：针对 $x$ 的梯度的停机准则，当当前步的梯度范数小于该值时认为该条件满足
opts.sigma ：增广拉格朗日系数
opts.gamma ： $x$ 更新的步长
opts.verbose ：不为 0 时输出每步迭代信息，否则不输出

if ~isfield(opts, 'maxit'); opts.maxit = 5000; end
if ~isfield(opts, 'sigma'); opts.sigma = 0.5; end
if ~isfield(opts, 'ftol'); opts.ftol = 1e-8; end
if ~isfield(opts, 'gtol'); opts.gtol = 1e-6; end
if ~isfield(opts, 'gamma'); opts.gamma = 1.618; end
if ~isfield(opts, 'verbose'); opts.verbose = 1; end

迭代准备。

tt = tic;
k = 0;
x = x0;
out = struct();

初始化对偶问题的对偶变量 $y$ 。

[m, ~] = size(A);
sm = opts.sigma;
y = zeros(m,1);

记录在初始时刻原问题目标函数值。

f = .5*norm(A*x - b,2)^2 + mu*norm(x,1);
fp = inf;
out.fvec = f;
nrmC = inf;

Cholesky 分解， $R$ 为上三角矩阵且 $R^\top R=A^\top A + \sigma I_m$ 。与原始问题同样的，由于罚因子在算法的迭代过程中未变化，事先缓存 Cholesky 分解可以加速迭代过程。

W = eye(m) + sm * (A * A');
R = chol(W);

迭代主循环

迭代主循环，当 (1) 达到最大迭代次数或 (2) 目标函数的变化小于阈值或 (3) 自变量 $x$ 的变化量小于阈值时，退出迭代循环。

while k < opts.maxit && abs(f - fp) > opts.ftol && nrmC > opts.gtol
    fp = f;

对 $z$ 的更新为向无穷范数球做欧式投影， $\displaystyle z^{k+1}=\mathcal{P}_{\|\cdot\|_\infty\le\mu} \left(x^k/\sigma-A^\top y^k\right)$ 。

    z = proj( - A' * y + x / sm, mu);

针对 $y$ 的子问题，即为求解线性方程组 $(I+\sigma AA^\top)y^{k+1}=A(x^k-\sigma z^{k+1})-b$ 。这里同样利用了事先缓存的 Cholesky 分解来加速 $y^{k+1}$ 的计算。

    h = A * (- z*sm + x) - b;
    y = R \ (R' \ h);

令 $c^{k+1}=A^\top y^{k+1} + z^{k+1}$ 为等式约束的约束违反度。增广拉格朗日函数对 $x$ 的梯度为 $\frac{\partial L_\rho(y,z,x)}{\partial x}=-\sigma c$ 。针对 $x$ 的子问题，进行一步梯度上升， $x^{k+1}=x^k - \gamma\sigma (A^\top y^{k+1} + z^{k+1})$ 。利用 nrmC （约束违反度的范数）作为停机判断依据。

    c = z + A' * y;
    x = x - opts.gamma * sm * c;
    nrmC = norm(c,2);

计算更新后的目标函数值，记录在 out.fvec 中。当 opts.verbose 不为 0 时输出详细的迭代信息。

    f = .5*norm(A*x - b,2)^2 + mu*norm(x,1);
    if opts.verbose
        fprintf('itr: %4d\tfval: %e\tfeasi: %.1e\n', k, f, nrmC);
    end
    k = k + 1;
    out.fvec = [out.fvec; f];
end

记录输出信息。

out.y = y;
out.fval = f;
out.itr = k;
out.tt = toc(tt);
out.nrmC = nrmC;

end

辅助函数

到无穷范数球 $\{x\big| \|x\|_\infty \le t\}$ 的投影函数。

function w = proj(x, t)
w = min(t, max(x, -t));
end

参考页面

在页面实例：交替方向乘子法解 LASSO 问题中我们展示此算法的一个应用。另外，对 LASSO 原问题的 ADMM 参考页面利用交替方向乘子法求解 LASSO 原问题。

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

版权声明