LASSOをまとめていると、LARS（最小角回帰法：Least Angle Regression）が出てきたので、これを機会に手順をちゃんと理解しようと思う。
行列と角度、高校や大学で習ったものの、いざ実データで、となるとどうしてもピンと来ない。実装してようやく理解できたので、その記録。
LARSとは
LARSは、LASSOを効率的に計算できるスパースモデリングの一種です。
回帰分析において、変数選択と回帰係数の推定を同時に行うアルゴリズムであり、特徴量が非常に多い高次元データセットに対して効率的に特徴選択できる強みがある。
LARS は最小二乗法をベースに次の流れで進みます。
- 初期値として全ての回帰係数を 0 にする。
- 目的変数となす角度が最も小さい（目的変数との相関係数の絶対値が最も大きい）説明変数を求め、この選択した説明変数をアクティブセットにする。
- アクティブセットの線形結合で表される方向に沿って、目的変数との残差が小さくなるように回帰係数を更新する。
- 更新した後の残差となす角度が最も小さい説明変数を求め、この選択した説明変数をアクティブセットに加える。
- アクティブセットの線形結合で表される方向（プロジェクションの方向）に沿って、残差を減らすように回帰係数を更新する。
- この残差とアクティブセット以外の列ベクトルとの角度が同じになるまで、4,5 を繰り返す。
- 残差とアクティブセット以外の列ベクトルとの角度が同じになったら、その列ベクトルをアクティブセットに加える。また、残差ベクトルとアクティブセット内の列ベクトルとの角度が 180 度になったら、その列ベクトルをアクティブセットから除く。
- 全ての説明変数がアクティブセットになるか、残差が 0 になったら終了する。
実行手順
この手順を書き下してみると、こんな流れになる。
lars <- function(X, y) { # 説明変数行列 `X` と目的変数ベクトル `y` が引数。 # 変数 `n` と `p` は、それぞれサンプル数と説明変数の数。 # `w` は回帰係数の初期値を保持する行列で、0が初期値。 # `mu` は予測値の初期値を保持するベクトルで、0 が初期値。 # `active_set` はアクティブな説明変数のインデックスを保持するベクトルで、空のベクトルが初期値。 # `residual` は残差を保持するベクトルで、初期値は目的変数 `y`。 n <- nrow(X) # サンプル数 p <- ncol(X) # 説明変数の数 w <- matrix(0, nrow = p, ncol = p) # 回帰係数の初期値 mu <- rep(0, n) # 予測値の初期値 active_set <- c() # アクティブな説明変数のインデックス residual <- y # 残差、初期値は目的変数 # 変数 `k` のループは、説明変数の数 `p` に対して実行されます。 for (k in 1:p) { # `k` が `p` より小さければ、アクティブな説明変数の係数と残差を更新 c <- t(X) %*% residual # 相関ベクトル c_max <- which.max(abs(c)) # 相関の絶対値が最大のインデックス active_set <- c(active_set, c_max) # `active_set` に `c_max` を追加 # アクティブな説明変数の部分行列を取得し、最小二乗法で回帰係数を推定 X_active <- X[, active_set, drop = FALSE] # アクティブな説明変数の部分行列 G <- solve(t(X_active) %*% X_active) %*% t(X_active) # Gは疑似逆行列 a <- drop(G %*% y) # 回帰係数の推定値 w[active_set, k] <- a # 係数行列に係数を追加（k列目に格納） # 予測値の更新 mu <- X_active %*% a # アクティブな説明変数のみを用いて予測値を計算 residual <- y - mu # 残差の更新 # アクティブセットのインデックス active_set、回帰係数 w 、予測値 mu、残差 residual を返す list(active_set = active_set, w = w, mu = mu, residual = residual) }
「最小二乗法」を振り返ったこと、「プロジェクションの方向」を幾何的に理解できたこと、この２つの理解がLARSの理解に大きくかかわりました。
