ぽんのブログ -22ページ目

ホームピグアメブロ

芸能人ブログ人気ブログ

ぽんのブログ

自分用の備忘録ブログです。書いてある内容、とくにソースは、後で自分で要点が分かるよう、かなり簡略化してます（というか、いい加減）。あまり信用しないように（汗

<< 前ページ次ページ >>

Lasso-Lars その８計算手順再掲

停止条件も含めて計算手順をまとめます。またそれぞれの操作を octave でも書いてみました。

ここで、前回 1A ( = ones (length(A), 1) ) としていた部分を s ( = sign (βj) = sign (cj) ) で置き換えます。それに伴い XA の部分も変更します。こちらの方がすっきりします。

観測データベクトルを y 、説明変数を列に格納した行列を X
X のサイズを m、 n に格納

[m, n] = size(X);

入力データ X、y を基準化及び中心化

y = y - mean(y);
for i = 1:n
Xi = X( : , i);
Xi= Xi - mean(Xi);
Xi = Xi / sqrt( Xi' * Xi);
end

１．初期値設定

$\boldsymbol{\mu}={\bf 0},\, \boldsymbol{\beta}={\bf 0},\, A=\phi$
停止フラグをリセット

mu = zeros(m, 1);    % μ = 0
beta = zeros(n , 1); % β= 0
A = [];                  % A = φ
stop_flag = false;

２．correlationの更新

${\bf c}={\bf X}^T\cdot({\bf y}- \boldsymbol{\mu})$

c = X' * (y - mu);

相関の絶対値の最大値を更新

$\hat{C}=\max_{j\in A}\{|c_j|\}$
A = φ ならその指標を記録
$i_{max}=i,\, i\in A,\, \hat{C}=|c_i|$

[ hat_c, c_amax_idx_tmp] = max( abs(c) );

if ( length(A) == 0 ),
c_amax_idx = c_amax_idx_tmp;
end

３．アクティブセットを更新

$i_{max}\geq 0$ ならそれをアクティブセットに加える。

$A=A+\{i_{max}\}$

if c_amax_idx >= 0,
A = [A, c_amax_idx];
end

４．等角ベクトルを更新

${\bf X}_A=(\cdots,{\bf x}_j,\cdots)_{j\in A}$

${\bf g}_A={\bf X}_A^T\cdot{\bf X}_A$

$s_j={\rm sign}(c_j),\, j\in A$

$A_A=1/\sqrt{{\bf s}^T\cdot{\bf g}_A^{-1}\cdot{\bf s}}$

${\bf w}_A=A_A{\bf g}_A^{-1}\cdot{\bf s}$

${\bf u}_A={\bf X}_A\cdot{\bf w}_A$

XA = X( : , A);                          % XA
GA = XA' * XA;                        % gA
invGA = inv(GA);
s = sign( c(A) );                     % s
absA = 1 / (s' * invGA * s);        % A_A
u = XA * wa;                           % uA

５．ステップサイズを更新

$\hat{\gamma}=\min_{j\in A_c}^+\left\{ \frac{\hat{C}-c_j}{A_A-a_j},\frac{\hat{C}+c_j}{A_A+a_j} \right\}$

ここに ${\bf a}={\bf X}^T\cdot{\bf u}_A$

a = X' * u;
Ac = setdiff(1 : n, A);     % complement of A
for i = 1:length(Ac),
j = Ac(i);
e1 = (hat_C - c(j) ) / (absA - a(j) );
e2 = (hat_C + c(j) ) / (absA + a(j) );
if e1 <= 0., e1 = Inf; end
if e2 <= 0., e2 = Inf; end
gamma_hat_tmp(i) = min(e1, e2);
end

[ gamma_hat, idx ] = min( gamma_hat_tmp );
gamma_hat_idx = Ac(idx);

$\widetilde{\gamma}=\min_{j\in A}^+\{-\beta_j/w_j\}$

gamma_tilde_tmp = - beta (A) . / w;
gamma_tilde_tmp( gamma_tilde_tmp <= 0 ) = Inf;
[ gamma_tilde, gamma_tilde_idx ] = min(gamma_tilde_tmp);
if isnan(gamma_tilde),
gamma_tilde = Inf;
end

$\hat{\gamma}<\widetilde{\gamma}$ なら $\gamma=\hat{\gamma}$
$\widetilde{\gamma}\leq \hat{\gamma}$ なら $\gamma=\widetilde{\gamma}$

if gamma_hat < gamma_tilde,
step_size = gamma_hat;
else
step_size = gamma_tilde;
end

６．βを更新

$\beta_j'=\beta_j+\gamma w_{Aj},\, j\in A$

beta_new = zeros(n, 1);
beta_new(A) = beta(A) + step_size * w;

７．βのL1ノルムと閾値 t の比較

$|\boldsymbol{\beta}|\leq t,\, t \leq|\boldsymbol{\beta}'|$

ならステップサイズを以下に更新

$\gamma=A_A(t-t_0),\, t_0=\sum_j|\beta_j|$
βを以下の線形補完値に更新しなおし、停止フラグを立てる。

$\beta_j'=\beta_j+\gamma w_{Aj},\, \, j\in A$

nrm1_prev = sum( abs( beta ) );
nrm1_cur = sum( abs( beta_new ) );

if nrm1_prev <= t && t <= nrm1_cur,
step_size = absA * (t - nrm1_prev);
beta_new(A) = beta(A) + step_size * w;
stop_flag = true;
end

８．β、μを更新

β = β'
$\boldsymbol{\mu}=\boldsymbol{\mu}+\gamma{\bf u}_A$

beta = beta_new;
mu = mu + step_size * u;

９．ストップフラグが立っていなければアクティブセットを更新

$\hat{\gamma}<\widetilde{\gamma}$ 　なら $\hat{\gamma}$ に対応する指標 j について $i_{max}=j$ とおく。
$\widetilde{\gamma}\leq \hat{\gamma}$ 　なら $\widetilde{\gamma}$ に対応する指標をアクティブセットから除く。

if ~stop_flag,
if gamma_hat < gamma_tilde,
    c_amax_idx = gamma_hat_idx;
else
    A( gamma_tilde_idx ) = [];
    c_amax_idx = -1;
end
end

以下 length(A) < n && !stop_flag の間２～９を繰り返す。

<< 前ページ次ページ >>