LASSO回帰を実装する
Rのmtcarsを使った実装例を作った。
回帰係数（w）と正則化パラメータ（$\lambda$）の初期値は適当に決めている。
回帰係数の初期値
ゼロベクトルや小さなランダムノイズを含んだベクトルなどがある。
- ゼロベクトル: 回帰係数をすべてゼロでスタートすると、最小二乗法の推定と同じくらいの大きさになる。
- ランダムノイズ: 小さなランダムノイズでスタートすると、初期値からバラバラになるので、局所最適解に陥る可能性が軽減する。
- データの特性に基づく初期値: データの特性などに基づいて選択する。例えば、特定の変数が目的変数と関連性が高いときにその係数を初期値にする、など。
正則化パラメータ
$\lambda$の大きさは、残る変数の数に影響する。$\lambda$が大きいと0になる係数が増え、変数は残らない。
$\lambda$の値は実際のデータやモデルの複雑さに依存する、つまり明確な設定基準はなさそう。
モデルの予測性能を探索して、その結果から$\lambda$を選択するようだ。
- クロスバリデーション（Cross-validation）: データを訓練データとテストデータに分割して、モデルを評価する。複数の$\lambda$に対して評価を繰り返して、予測性能が最も高い$\lambda$を選ぶ。
- 情報量基準（Information Criterion）: AIC（赤池情報量基準）やBIC（ベイズ情報量基準）などを使用して、モデルの適合度とモデルの複雑さのバランスが良い$\lambda$を選ぶ。
- パス探索（Path Searching）: 一連の$\lambda$の値でのでのLasso回帰の推定結果を求めて、変数の選択状況や回帰係数の変化を視覚的に確認して、適切な$\lambda$を選ぶ。
回帰係数wの求め方
i番目の回帰係数は、説明変数と残差の内積と$\lambda$の大きさできまる。
この内積が$\pm \frac{\lambda}{2}$に入れば、$\lambda$と説明変数の大きさで調整した値、入らなければ0になる。
もちろん、説明変数が全て0なら、回帰係数は0になる。
説明変数X[,i]と残差rの内積 |
回帰係数w[i] |
|---|---|
| 説明変数と残差の内積 < $\frac{\lambda}{2}$ | $\lambda$ の半分を足して説明変数の二乗和で割る |
| 説明変数と残差の内積 < -$\frac{\lambda}{2}$ | $\lambda$ の半分を引いて説明変数の二乗和で割る |
| その外側 | 0 |
| 説明変数が全て0 | 0 |
lasso <- function(X, y, w, lambda) { # 座標降下法（Coordinate Descent）でLasso回帰の係数を計算する関数 ## X: 説明変数のマトリクス ## y: 目的変数のベクトル ## lambda: L1正則化のλ # ラムダの半分 lambda2 <- lambda / 2 for (i in 1:ncol(X)) { # 残差の二乗和とL1正則化項の和をrにする r <- sum((y - X %*% w)^2) + lambda * sum(abs(w)) if (is.na(r) | is.infinite(r)) r <- 0 # rを計算できないときのエラー処理 # 説明変数の列と残差の内積 sumXr <- sum(X[, i] * r) # 説明変数の列の二乗和 sumX2 <- sum(X[, i]^2) if (sumX2 == 0) { w[j] <- 0 } else { if (sumXr > lambda2) { w[i] <- (sumXr - lambda2) / sumX2 } else if (sumXr < -lambda2) { w[i] <- (sumXr + lambda2) / sumX2 } else { w[i] <- 0 } } } return(w) } data(mtcars) X <- as.matrix(mtcars[, -1]) y <- mtcars$mpg w = rep(0, ncol(X)) names(w) <- colnames(X) lambda <- .0000000001 lasso(apply(X,2,scale), y, w, lambda)
これで、L1正則化を利用した回帰係数が分かる。
でも実際は、目的変数と予測値の残差が小さくなる（収束する）まで繰り返すことになる。
実際は、最大繰り返し数（max_iter）と収束値（tol）を組み込むことになる。
lasso_iter <- function(X, y, lambda, max_iter = 1000, tol = lambda/1000) { # 座標降下法（Coordinate Descent）でLasso回帰の係数を計算する関数 ## X: 説明変数のマトリクス ## y: 目的変数のベクトル ## lambda: L1正則化のλ # 収束判定するための引数 ## max_iter: 繰り返し回数の上限 ## tol: 収束判定の下限値 # 回帰係数の初期値 w <- rep(0, ncol(X)) names(w) = colnames(X) # 収束の確認フラグ converged <- FALSE # 反復回数のカウント iter <- 0 # 収束するか最大反復回数に達するまで繰り返す while (!converged && iter < max_iter) { # 前の回帰係数のベクトルを保存する w_old <- w # 各説明変数に対して回帰係数を更新する w = lasso(X, y, w, lambda) # 反復回数を増やす iter <- iter + 1 # 前の回帰係数と現在の回帰係数の差分ベクトルのノルムを計算する diff_norm <- sqrt(sum((w - w_old)^2)) # 差分ベクトルのノルムが許容誤差より小さくなったら収束したと判断する if (diff_norm < tol) { converged <- TRUE } } # 収束判定、反復回数、回帰係数を返す return(list(converged = converged, iter = iter, w = w)) } data(mtcars) X <- as.matrix(mtcars[, -1]) y <- mtcars$mpg w = rep(0, ncol(X)) names(w) <- colnames(X) lambda <- .0000000001 lasso_iter(apply(X, 2, scale), y, lambda)
