ADMM解LASSO問題
本文為《最優化：建模、算法與理論》8.6 交替方向乘子法 對應的代碼及其說明代碼作者：文再文、劉浩洋、戶將為方便閱讀將代碼注釋中
本文為《最優化：建模、算法與理論》8.6 交替方向乘子法 對應的代碼及其說明
代碼作者：文再文、劉浩洋、戶將
為方便閱讀將代碼注釋中的latex寫出來了
利用 ADMM 求解 LASSO 問題的原問題
針對 LASSO 問題
引入拉格朗日乘子 
在 ADMM 的每一步迭代中，交替更新 
%% 初始化和迭代準備
% 函數通過優化上面給出的增廣拉格朗日函數，以得到 LASSO 問題的解.
% 輸入信息： 
% 輸出信息： 迭代得到的解 
% * |out.fvec| ：每一步迭代的 LASSO 問題目標函數值
% * |out.fval| ：迭代終止時的 LASSO 問題目標函數值
% * |out.tt| ：運行時間
% * |out.itr| ：迭代次數
% * |out.y| ：迭代終止時的對偶變量
% * |out.nrmC| ：約束違反度，在一定程度上反映收斂性
function [x, out] = LASSO_admm_primal(x0, A, b, mu, opts)
%%%
% 從輸入的結構體 |opts| 中讀取參數或采取默認參數.
% * |opts.maxit| ：最大迭代次數
% * |opts.ftol| ：針對函數值的停機準則，當相鄰兩次迭代函數值之差小于該值時認為該條件滿足
% * |opts.gtol| ：針對
% * |opts.sigma| ：增廣拉格朗日系數
% * |opts.gamma| ：
% * |opts.verbose| ：不為 0 時輸出每步迭代信息，否則不輸出
if ~isfield(opts, 'maxit'); opts.maxit = 5000; endnif ~isfield(opts, 'sigma'); opts.sigma = 0.01; endnif ~isfield(opts, 'ftol'); opts.ftol = 1e-8; endnif ~isfield(opts, 'gtol'); opts.gtol = 1e-14; endnif ~isfield(opts, 'gamma'); opts.gamma = 1.618; endnif ~isfield(opts, 'verbose'); opts.verbose = 1; end
%%%
% 迭代準備
k = 0;ntt = tic;nx = x0;nout = struct();
%%%
% 初始化 ADMM 的輔助變量 
[m, n] = size(A);nsm = opts.sigma;ny = zeros(n,1);nz = zeros(n,1);
%%%
% 計算并記錄起始點的目標函數值。
fp = inf; nrmC = inf;nf = Func(A, b, mu, x);nf0 = f;nout.fvec = f0;
%%%
% Cholesky 分解， 
AtA = A'*A;nR = chol(AtA + opts.sigma*eye(n));nAtb = A'*b;
%% 迭代主循環
% 迭代主循環，當 (1) 達到最大迭代次數或 (2) 目標函數的變化小于閾值或 (3) 自變量
while k < opts.maxit && abs(f - fp) > opts.ftol && nrmC > opts.gtoln fp = f;
%%%
% 更新 
% 
w = Atb + sm*z - y;n x = R (R' w);
%%%
% 更新 
% 即 
c = x + y/sm;n z = prox(c,mu/sm);
%%%
% 以 
% 更新 
% 以 
y = y + opts.gamma * sm * (x - z);n f = Func(A, b, mu, x);n nrmC = norm(x-z,2);
%%%
% 輸出每步迭代的信息. 迭代步 
if opts.verbosen fprintf('itr: %4dtfval: %etfeasi:%.1en', k, f,nrmC);n endn k = k + 1;n out.fvec = [out.fvec; f];nend
%%%
% 退出循環，記錄輸出信息。
out.y = y;nout.fval = f;nout.itr = k;nout.tt = toc(tt);nout.nrmC = norm(c - y, inf);nend
%% 輔助函數
%%%
% 函數 
function y = prox(x, mu)ny = max(abs(x) - mu, 0);ny = sign(x) .* y;nend
%%%
% LASSO 問題的目標函數 
function f = Func(A, b, mu, x)nw = A * x - b;nf = 0.5 * (w' * w) + mu*norm(x, 1);nend
利用 ADMM 求解 LASSO 對偶問題
針對 LASSO 問題
考慮其對偶問題的 ADMM 等價形式
引入
在 ADMM 的每一步迭代中，交替更新 
%% 初始化和迭代準備
% 函數通過優化上面給出的增廣拉格朗日函數，以得到 LASSO 問題的解。
% 輸入信息：
% 輸出信息： 迭代得到的解
% * |out.fvec| ：每一步迭代的 LASSO 問題目標函數值
% * |out.fval| ：迭代終止時的 LASSO 問題目標函數值
% * |out.tt| ：運行時間
% * |out.itr| ：迭代次數
% * |out.y| ：迭代終止時的對偶變量
% * |out.nrmC| ：迭代終止時的約束違反度，在一定程度上反映收斂性
function [x, out] = LASSO_admm_dual(x0, A, b, mu, opts)
%%%
% 從輸入的結構體 |opts| 中讀取參數或采取默認參數。
% * |opts.maxit| ：最大迭代次數
% * |opts.ftol| ：針對函數值的停機準則，當相鄰兩次迭代函數值之差小于該值時認為該條件滿足
% * |opts.gtol| ：針對
% * |opts.sigma| ：增廣拉格朗日系數
% * |opts.gamma| ：
% * |opts.verbose| ：不為 0 時輸出每步迭代信息，否則不輸出
if ~isfield(opts, 'maxit'); opts.maxit = 5000; endnif ~isfield(opts, 'sigma'); opts.sigma = 0.5; endnif ~isfield(opts, 'ftol'); opts.ftol = 1e-8; endnif ~isfield(opts, 'gtol'); opts.gtol = 1e-6; endnif ~isfield(opts, 'gamma'); opts.gamma = 1.618; endnif ~isfield(opts, 'verbose'); opts.verbose = 1; end
%%%
% 迭代準備。
tt = tic;nk = 0;nx = x0;nout = struct();
%%%
% 初始化對偶問題的對偶變量
[m, ~] = size(A);nsm = opts.sigma;ny = zeros(m,1);
%%%
% 記錄在初始時刻原問題目標函數值。
f = .5*norm(A*x - b,2)^2 + mu*norm(x,1);nfp = inf;nout.fvec = f;nnrmC = inf;
%%%
% Cholesky 分解， 
% 與原始問題同樣的，由于罰因子在算法的迭代過程中未變化，事先緩存 Cholesky 分解可以加速迭代過程。
W = eye(m) + sm * (A * A');nR = chol(W);
%% 迭代主循環
% 迭代主循環，當 (1) 達到最大迭代次數或 (2) 目標函數的變化小于閾值或
% (3) 自變量
while k < opts.maxit && abs(f - fp) > opts.ftol && nrmC > opts.gtoln fp = f;
%%%
% 對 
z = proj( - A' * y + x / sm, mu);
%%%
% 針對 
% 這里同樣利用了事先緩存的 Cholesky 分解來加速 
h = A * (- z*sm + x) - b;n y = R (R' h);
%%%
% 令 
% 增廣拉格朗日函數對 
% 針對 
% 利用 |nrmC| （約束違反度的范數）作為停機判斷依據。
c = z + A' * y;n x = x - opts.gamma * sm * c;n nrmC = norm(c,2);
%%%
% 計算更新后的目標函數值，記錄在 |out.fvec| 中。當 |opts.verbose| 不為 0 時輸出詳細的迭代信息。
f = .5*norm(A*x - b,2)^2 + mu*norm(x,1);n if opts.verbosen fprintf('itr: %4dtfval: %etfeasi: %.1en', k, f, nrmC);n endn k = k + 1;n out.fvec = [out.fvec; f];nend
%%%
% 記錄輸出信息。
out.y = y;nout.fval = f;nout.itr = k;nout.tt = toc(tt);nout.nrmC = nrmC;nend
%% 輔助函數
% 到無窮范數球 
function w = proj(x, t)nw = min(t, max(x, -t));nend
考慮LASSO問題：
首先考慮利用 ADMM 求解原問題：將其轉化為 ADMM 標準問題
則可以利用 ADMM 求解. 相應的，對于 LASSO 對偶問題
則等價于
對于上述的兩個等價問題利用 ADMM 求解.
%% 構造 LASSO 問題
% 設定隨機種子.
clear;nseed = 97006855;nss = RandStream('mt19937ar','Seed',seed);nRandStream.setGlobalStream(ss);
%%%
% 構造 LASSO 優化問題
% 
% 生成隨機的矩陣 
m = 512;nn = 1024;nA = randn(m, n);nu = sprandn(n, 1, 0.1);nb = A * u;nx0 = randn(n, 1);nmu = 1e-3;
%% 利用 ADMM 求解 LASSO 問題
% 首先在更嚴格的停機準則下進行試驗，將收斂時得到的歷史最優函數值作為真實的最優值的參考 
opts = struct();nopts.verbose = 0;nopts.maxit = 2000;nopts.sigma = 1e-2;nopts.ftol = 1e-12; nopts.gtol = 1e-15;n[x, out] = LASSO_admm_primal(x0, A, b, mu, opts);nf_star = min(out.fvec);
%%%
% 利用 ADMM 求解 LASSO 對偶問題.
opts = struct();nopts.verbose = 0;nopts.maxit = 2000;nopts.sigma = 1e2;nopts.ftol = 1e-8; nopts.gtol = 1e-10;n[x, out] = LASSO_admm_dual(x0, A, b, mu, opts);ndata1 = (out.fvec - f_star)/f_star;nk1 = length(data1);
%%%
% 利用 ADMM 求解 LASSO 原問題.
opts = struct();nopts.verbose = 0;nopts.maxit = 2000;nopts.sigma = 1e-2;nopts.ftol = 1e-8; nopts.gtol = 1e-10;n[x, out] = LASSO_admm_primal(x0, A, b, mu, opts);ndata2 = (out.fvec - f_star)/f_star;nk2 = length(data2);
%% 結果可視化
% 對每一步的目標函數值與最優函數值的相對誤差進行可視化.
fig = figure;nsemilogy(0:k1-1, data1, '-', 'Color',[0.2 0.1 0.99], 'LineWidth',2);nhold onnsemilogy(0:k2-1, data2, '-.','Color',[0.99 0.1 0.2], 'LineWidth',1.5);nlegend('ADMM解對偶問題','ADMM解原問題');nylabel('$(f(x^k) - f^*)/f^*$', 'fontsize', 14, 'interpreter', 'latex');nxlabel('迭代步');nprint(fig, '-depsc','admm.eps');
結果分析
上圖反映了使用 ADMM 求解 LASSO 原問題和對偶問題的表現。在兩個問題上目標函數值都存在波動，表明 ADMM 并非下降算法，然而在沒有使用連續化策略的情況下，ADMM 依然達到了收斂，這說明 ADMM 較之其它算法具有一定的優越性。注意，雖然在這一例子中 ADMM求解原問題所需迭代次數較少，但求解對偶問題時每一步迭代時間更短。綜合看來在該例子中 ADMM 求解對偶問題的性能更高。
