Coordinate descent for lasso ３.CovarianceUpdates

Lasso の場合の coordinate descent での解の更新式は、前回載せたように



でした。S の左の項は更に以下の様に書きかえる事が出来ます。



つまり

式（４）


となります。ここに で、これは計算中不変です。
式（４）の方が計算量は少なくなりますね。

と言う訳で基本的な考え方は以上で、とてもシンプルです。
実装も（細かい事を考えなければ）とても簡単です。
以下は elastic net を (cyclic) coordinate descent で解くスクリプトです。
diabetes.data のソリューションパスを求めています。

=== soft_threshold.m ===
function val = soft_threshold (z, gamma)
    val = 0.;
    if gamma < abs(z),
        if z > 0,
            val = z - gamma;
        else
            val = z + gamma;
        end
    end
end

続いて coordinate descent のカーネル部分

=== cdescent.m ===

function [beta, mu] = cdescent_cyclic(y,X,c,beta,mu,lambda1,lambda2,maxiter)

    [n, p] = size(X);
    tol = 1.e-5;

    if isempty(beta),
        beta = zeros(p, 1);
        mu = zeros(n, 1);
    end

    if isempty(c),
        c = X' * y;
    end

    iter = 1;
    while iter < maxiter,
        delta = beta; % beta の差分: delta = beta_prev で初期化
        for j = 1:p,
            beta(j, 1) = soft_threshold(
　　　　　　　　　　　　　　c(j) - X(:, j)' * mu + beta(j),
　　　　　　　　　　　　　　lambda1) / (1 + lambda2);
 
           delta(j) = delta(j) - beta(j); % = beta_prev(j) - beta(j)
            if abs(delta(j)) > 0, % mu を更新
                mu = mu - X(:, j) * delta(j);
            end
        end

        if sum(abs(delta)) < tol, % 収束判定
            break;
        end

        iter = iter + 1;
    end
end

j = 1 : p のループの中で、soft thresholding で各々の β_j を更新しています。
β_j が更新されたら μ = X β もそれに合わせ更新します。
更新後の解 β' は β_j だけが



と更新されたとします。すると



となります。

で、全ての β を更新し、前のβから変化していなければ収束したものとして while ループを抜ける、そうでなければ j = 1 : p のループをやり直す・・・とします。

以上の関数で diabetes.data を解析するのか下のスクリプトです。

=== test_cdescent.m ===

clear all;

data = dlmread("diabetes.data", '\t', 1, 0);
y = data(:, end);
X = data(:, 1:end-1);

y = y - mean(y);
X = bsxfun(@minus, X, mean(X));
X = bsxfun(@rdivide, X, sqrt(dot(X, X)));

lambda2 = 1;

[n, p] = size(X);

b = zeros(p, 1);
mu = zeros(n, 1);
c = X' * y;
lmax = max(abs(c));    % lambda1_max
imax = round(log10(lmax));

t = [];
beta = [];
for i = imax:-0.1:imax-3
    lambda1 = 10^i;
    [b, mu] = cdescent_cyclic(y, X, c, b, mu, lambda1, lambda2, 1000);
    t = [t, lambda1];
    beta = [beta, b];
end

% elastic net solution
beta = beta * (1 + lambda2);
plot (sum(abs(beta)), beta, '-x')

このスクリプトを実行し |β| を横軸にして書いたソリューションパスがこちらです。