Lasso回帰の最適化アルゴリズム
最適化アルゴリズムは座標降下法が有名ですが、調べてみると他にもいくつかありました。 分かる範囲でリストアップしてみた。 それぞれの理解をまとめておく。
| 名称 | 特徴 | メリット | デメリット |
|---|---|---|---|
| 座標降下法 | 1変数ずつ重みの更新を繰り返して 最適解に近似させる |
計算が簡単で高速 | 収束性が保証されない |
| サイクリック座標降下法 | 座標降下法の一種 重みを固定した順番で更新する |
座標降下法よりも収束が速い | 更新順序によっては収束しない |
| シューティングアルゴリズム | 座標降下法の一種 重みをランダムに選んで更新する |
重みの更新順序に依存しない | 収束までに多くの反復が必要 |
| LARS（最小角回帰法） | 重みを最も急速に増加させる 方向に沿って更新する |
変数選択と正則化を同時に行える | 計算量が大きい |
| プロキシミティ演算子法 | 重みを現在の値から一定量だけ動かす | 非線形回帰や非凸最適化問題にも適用できる | 収束速度が遅い |
関数にする
数式を見ていても、今ひとつピンとこないので、動作の理解も含めてコーディングしてみた。
座標降下法（Coordinate Descent）
- 残差の二乗和とL1正則化項を計算し、
rとして保存する。 - もし説明変数の列がすべて0の場合は、回帰係数を0にする。
- そうでなければ、回帰係数を以下の条件で計算する。
ゼロを中心にした$\lambda$の範囲外なら回帰係数はゼロ、範囲内なら非ゼロという感じです。
説明変数X[,i]と残差rの内積 |
回帰係数w[i] |
|---|---|
| 説明変数と残差の内積 < $\lambda$の半分 | $\lambda$ の半分を足して説明変数の二乗和で割る |
| 説明変数と残差の内積 < -$\lambda$の半分 | $\lambda$ の半分を引いて説明変数の二乗和で割る |
| その外側 | 0 |
| 説明変数が全て0 | 0 |
lasso <- function(X, y, w, lambda) { # ラムダの半分 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) }
サイクリック座標降下法（Cyclic Coordinate Descent）
- 残差を
rにする。 - 説明変数の列ごとに次の手順を繰り返す。
- 「説明変数の列と残差の内積」を
zにする。 - 既存の回帰係数
w[-i]の符号とlambda/nrow(X)を掛けた値をaにする。 - もし、zの絶対値が
lambda/nrow(X)より大きければ、回帰係数をaにする。 - それ以外なら、回帰係数を
0にする。
- 「説明変数の列と残差の内積」を
lasso <- function(X, y, w, lambda) { for (i in 1:ncol(X)) { r <- y - X %*% w z <- X[, i] %*% r / nrow(X) a <- z + sign(w[-i]) * lambda / nrow(X) w[j] <- ifelse(abs(z) > lambda / nrow(X), a, 0) } return(w) }
シューティングアルゴリズム（Shooting Algorithm）
- 「説明変数の列と目的変数の内積」の符号から回帰係数の符号を決める。
abs(x * y) - lambda / 2を0以下に制約した結果を、max(abs(x * y) - lambda / 2, 0)にする。- この値を説明変数の列の二乗和
sum(x^2)で割り、回帰係数にする。
lasso <- function(X, y, w, lambda) { for (i in 1:ncol(X)) { pm <- sign(sum(X[, i] * y)) w[i] <- pm * max(abs(X[, i] * y) - lambda / 2, 0) / sum(X[, i]^2) } return(w) }
最小角回帰法（LARS: Least Angle Regression）
- 最も相関の高い説明変数を見つける。
- アクティブな説明変数の部分行列を取得し、最小二乗法で回帰係数を計算する。
- アクティブな説明変数の係数と残差を更新する。
- アクティブな説明変数の係数のノルムを計算し、プロジェクション方向を決定する。
LARSは理解があいまいなので、別でまとめる。
lars <- function(X, y) { n <- nrow(X) # サンプル数 p <- ncol(X) # 説明変数の数 w <- matrix(0, nrow = p, ncol = p) # 回帰係数の初期値 mu <- rep(0, n) # 予測値の初期値 active_set <- c() # アクティブな説明変数のインデックス residual <- y # 残差、初期値は目的変数 for (k in 1:p) { # 最も相関の高い説明変数を見つける c <- t(X) %*% residual # 相関ベクトル c_max <- which.max(abs(c)) # 相関の絶対値が最大のインデックス active_set <- c(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) }
プロキシミティ演算子法
- 入力ベクトル
vと正則化項の重みパラメータlambdaを受け取る。 - L1正則化に基づくプロキシミティ演算子を計算する。
- 要素の絶対値から正則化項の重みを引く。
- 入力ベクトルの要素ごとに、ゼロ以下にする。
lasso <- function(v, lambda) { sign(v) * pmax(abs(v) - lambda, 0) }
