genlasso

パッケージの概要

genlassoパッケージはラッソ正則化線形回帰モデルについて、 通常のものよりも罰則項を一般化したモデルを取り扱うものです。 加えて、フューズドラッソ（Fused LASSO）等の特別なパターンについてはより効率的なアルゴリズムを用意しています。 なお、GLM には対応していません。

一般化ラッソ回帰問題とは

まず、本パッケージが取り扱う一般化ラッソ回帰問題について簡単に説明します。

実数値の説明変数が k 個・目的変数が1個あり、 これらについて n 組の観測がある場合の線形回帰問題を考えます。

通常の線形回帰問題は、説明変数の実測値を表す n 行 k 列実数値行列 X = (x_{ij}) と目的変数の実測値を表す y = (y_i)_{i=1}^n \in \mathbb{R} ^ n に対して、

\sum_{i=1}^{n}\left(y_i - \sum_{j=1}^{k} \left( x_{ij} \cdot \beta_j \right) \right)^2 を最小化するような \beta = (\beta_j)_{j=1}^k \in \mathbb{R}^k を求める問題（最小二乗法）として定式化されます。

ベクトル v = (v_i)_{i=1}^m \in \mathbb{R}^m と 1 \leq p < \infty に対し、 p - ノルム \|v\|_p := \left( \sum_{i=1}^m |v_i|^p \right)^{1/p} が定義されます。これを用いると、最小化すべき関数は

\| y - X\beta \|_2^2

と簡単に書けます。

これを少し変更した、

\frac12\| y - X\beta \|_2^2 + \lambda\| \beta \|_1

を最小化する問題をラッソ回帰問題 ¹ といいます。 ここで、\lambda \geq 0 はモデルのハイパーパラメータです。

この第2項をラッソ罰則項といい、これにより線形回帰モデルの解が安定しやすくなります。 このように、最小化問題にラッソ罰則項を付け加えることをラッソ正則化といい、 特に予測モデル構築の観点では、予測精度が向上し、かつ変数選択が自動的に行われることが知られています(岩沢 宏和 and 平松 雄司 2019)。 入力データに比して説明変数が非常に多い場合のような、いわゆるスパースモデリングのための重要手法として深く研究されています。

本パッケージで取り扱う一般化ラッソ回帰問題とは、この罰則項をさらに一般化した

\frac12\| y - X\beta \|_2^2 + \lambda\| D \beta \|_1 の最小化問題を取り扱うものです。ただし、D は l 行 k 列実数値行列です。

この中には、罰則項を \lambda \sum_{i=1}^{k-1} | \beta_{i+1} - \beta_{i} | とした、 いわゆるフューズドラッソ回帰問題も含まれます。 隣接した説明変数同士の係数が同じ値に近づくような作用を与え、ノイズ除去などに用いられます。 一般化ラッソ回帰問題にて D を次のように設定することでこの問題に帰着します。

D = \begin{bmatrix} -1 & 1 & 0 & 0 & \cdots & 0 & 0\\ 0 & -1 & 1 & 0 & \cdots & 0 & 0 \\ & & \vdots & & \ddots & \vdots & \\ 0 & 0 & 0 & 0 & \cdots & -1 & 1 \\ \end{bmatrix}

本パッケージが優れているのは、この一般化ラッソ回帰問題について Taylor B. Arnold and Ryan J. Tibshirani (2016) の手法を実装することにより、 解パス、すなわちすべての \lambda に対する解を一度に与えることができることです。 これにより、\lambda を変更したときの振る舞いを観察したり、 ハイパーパラメータ \lambda のチューニングに活用したりすることができます。

準備

パッケージの読み込み

library(genlasso)
library(ggplot2) #可視化に使用

データセットの読み込み

Bostonデータセットとはボストンの住宅価格に関するデータで、 今回はmedv（住宅価格）を目的変数とし、これを他の説明変数で予測するモデルを構築することとします。

set.seed(42)
#学習用データと評価用データを分割します。
df_all <- MASS::Boston
split_df_all <- rsample::initial_split(df_all)

#学習用データを説明変数と目的変数に分割します。
df_train <- rsample::training(split_df_all)
df_train_X <- dplyr::select(df_train, -medv)
df_train_y <- df_train$medv

#評価用データを説明変数と目的変数に分割します。
df_test <- rsample::testing(split_df_all)
df_test_X <- dplyr::select(df_test, -medv)
df_test_y <- df_test$medv

使用方法

通常のラッソ回帰問題

Bostonデータセットにて、罰則項を \lambda \|\beta\| とする通常のラッソ回帰モデルを構築することを考えます。

罰則項の存在から、説明変数のスケールの違いにより最小化問題の解が変わるという性質があります。

例えば、ある説明変数だけを全ての観測で一様に10倍すると、 通常の線形回帰問題では対応する係数を10分の1すればよいだけですが、 ラッソ回帰問題では罰則項における評価が他の説明変数に比べて10分の1になってしまい、 最小化問題の解が変わってしまいます。

スケールの違いによって解が変わらないようにするため、 平均がゼロ、分散を1になるように線形変換する標準化（normalize）を最初に行うこととします。

#標準化のため、recipesパッケージを使用します。
#「標準化する」という手順書を作成し、これを変数に格納します。
rec_nm <- recipes::recipe(df_train, medv ~ .) |> 
  #全ての説明変数を標準化
  recipes::step_normalize(recipes::all_numeric_predictors())

#学習用データに基づき、標準化に使用する変換方法を具体的に決定します。
#（これにより、評価用データに対しても学習データと同じルールで変換できます）
recp_nm <- rec_nm |> recipes::prep()

#標準化したデータを別の変数に格納します。
df_train_nm <- recp_nm |> recipes::bake(new_data = df_train)
df_train_X_nm <- dplyr::select(df_train_nm, -medv)

df_test_nm <- recp_nm |> recipes::bake(new_data = df_test)
df_test_X_nm <- dplyr::select(df_test_nm, -medv)

また、genlassoパッケージには切片項を明示的に追加する機能がないため、 学習用データに全件1を入力した列を追加することで対応します。

model.matrix関数は計画行列を作成するための汎用的な関数ですが、 今回の場合は以下のようにすることで、切片項の列を追加するのに使用できます。

df_train_1X_nm <- model.matrix(~ ., df_train_X_nm)

次に、罰則項の係数行列 D の設定をします。 通常のラッソ回帰の場合は単位行列 I を指定した場合に対応しますが、 切片項の分だけ D の列数も1列増やしておく必要があります。

D <- diag(1, ncol(df_train_X)+1)

さらに、切片項に対しては罰則項の対象外とするため、1行目は取り除きます。

D <- D[-1,]

これで、回帰問題の設定に用いる変数は全て用意できました。

後は、genlasso関数に目的変数の実測値y、説明変数の実測値（計画行列）X、 罰則項の係数行列Dを入力することで回帰モデルが構築できます。

model_genlasso <- genlasso(y = df_train_y,
                           X = df_train_1X_nm, 
                           D = D)

構築したモデルに対してplot関数を適用すると、 λによってどのように係数が変化するか（解パス）をプロットすることができます。

plot(model_genlasso)

このプロットにより、正則化によってどのように変数選択されていくかを可視化することができます。

今回のようにモデリングで使用する部分がグラフの左の方に固まってしまうケースでは その部分を拡大して観察したほうが良いように思われますが、 標準ではこのグラフを拡大する方法は用意されていません。

グラフの範囲は $lambda[1] によって決定されているので、 これを上書きすることで強引に拡大することは可能です。

model_genlasso_tmp <- model_genlasso
model_genlasso_tmp$lambda[1] <- 1*nrow(df_train_X)
plot(model_genlasso_tmp)

model_genlasso_tmp$lambda[1] <- 0.1*nrow(df_train_X)
plot(model_genlasso_tmp)

また、λと説明変数の個数の関係はsummary関数でも確認できます。

summary(model_genlasso)

   df    lambda     rss
    1   2511.89   30958
    2   2082.40   25737
    3   1072.16   15226
    4    404.69   10990
    5    294.27   10549
    6    291.33   10538
    7    248.56   10355
    8    173.02    9839
    9     94.48    9210
   10     94.16    9208
   11     87.99    9131
   12     54.55    8764
   13     30.97    8507

特定のλでの係数を確認するには、coef関数を使用します。

coef(model_genlasso, lambda = c(0.01, 0.1, 1, 10)*nrow(df_train_X))

$beta
             3.79          37.9           379          3790
 [1,] 22.44828496  2.244828e+01  2.244828e+01  2.244828e+01
 [2,] -1.07583207 -7.910061e-01  9.508378e-16  4.149459e-15
 [3,]  0.82295065  4.282012e-01  3.012037e-16 -6.106227e-15
 [4,] -0.03736614 -1.739666e-01  3.363591e-16  3.171075e-15
 [5,]  0.47094307  4.555891e-01 -7.536400e-16  2.844947e-16
 [6,] -2.10978311 -1.679800e+00 -1.552797e-15  4.274359e-15
 [7,]  2.49902980  2.618992e+00  2.546729e+00 -1.122713e-14
 [8,] -0.16356541 -5.582320e-16  2.360366e-15  5.152129e-15
 [9,] -3.24170209 -2.466214e+00  2.635306e-16  4.676814e-15
[10,]  2.63338674  1.281014e+00 -3.570671e-15  1.942890e-15
[11,] -1.74110608 -5.748326e-01  3.045836e-15 -7.410739e-15
[12,] -2.15187674 -2.025676e+00 -1.220700e+00 -5.925815e-15
[13,]  0.68719108  6.315578e-01  5.366666e-02 -2.997602e-15
[14,] -3.84759123 -3.899315e+00 -3.662125e+00 -2.975051e-15

$lambda
[1]    3.79   37.90  379.00 3790.00

$df
[1] 14 13  5  1

この回帰モデルで予測を行うには、predict関数を使用します。

#λの値はそれなりの精度になるものを手で検証して設定しました。
lambda_gen <- 0.027 * nrow(df_train_X)
preds <- predict(
  object = model_genlasso,
  Xnew = model.matrix(~ ., df_test_X_nm),
  lambda = lambda_gen #λの値を引数lambdaに与えます。
)
preds <- c(preds$fit)# 予測値を取り出します。

#精度評価のため、RMSEを計算します。
.rmse <- round(yardstick::rmse_vec(df_test_y, preds), 3)

#横軸に目的変数の実測値、縦軸に予測値をプロットします。
ggplot(mapping = aes(x = df_test_y, y = preds)) +
  geom_point() +
  ggtitle(paste("Lasso @", lambda_gen, " | RMSE:", .rmse))

λのチューニング

予測モデルとして用いる場合、 ラッソ回帰における \lambda を決定する方法として交差検証法（CV:Cross Validation）が考えられます。 今回はこれを実装してみましょう。

本パッケージの特徴として、1回のモデル構築で全てのλに対する解が得られるため、 λの候補数だけモデル構築を繰り返す必要はありません。 例えば 5-fold 交差検証法の場合は、累計5回モデル構築を行えば十分です。

set.seed(2024)
#k-fold 交差検証法に用いるデータセットを用意します。
#学習用データを5等分したうえで、
#そのうち4個（分析セット）でモデル構築 → 1個（検証セット）で精度確認　を5回繰り返すものですが、
#以下の関数により、分析セット・検証セットの組を計5組用意することができます。
splits_cv <- rsample::vfold_cv(df_train, v = 5)

#各foldにおける予測モデルを格納するリスト
models_cv <- list()
#各fold, λにおける結果を格納するデータフレーム
results_cv <- tibble::tibble(id = character(), lambda = numeric(), rmse = numeric())
#λの候補を設定します。
lambdas <- c(0, 10^(seq(-100, 0, 1)*3/100)) * nrow(df_train)

#5回モデル構築を繰り返します。
for(i in 1:nrow(splits_cv)){
  #i番目の分析セットを抽出します。
  split <- splits_cv$splits[[i]]
  df_analysis <- rsample::analysis(split)
  df_analysis_nm <- recp_nm |> recipes::bake(new_data = df_analysis)
  df_analysis_X_nm <- dplyr::select(df_analysis_nm, -medv)
  df_analysis_y <- df_analysis$medv

  #i版目の分析セットでモデルを構築し、出来上がったモデルをリストに格納します。
  models_cv[[splits_cv$id[[i]]]] <- 
    genlasso(y = df_analysis_y,
             X = model.matrix(~ ., df_analysis_X_nm), 
             D = D)
}

構築した5つのモデルに対して、事前に用意したλの候補に対する精度評価を行います。 評価指標（ここではRMSEを使用します）の平均が良好なλを、最終的なλとして採用します。

for(i in 1:nrow(splits_cv)){
  #i番目の検証セットを抽出します。
  split <- splits_cv$splits[[i]]
  df_assessment <- rsample::assessment(split)
  df_assessment_nm <- recp_nm |> recipes::bake(new_data = df_assessment)
  df_assessment_X_nm <- dplyr::select(df_assessment_nm, -medv)
  df_assessment_1X_nm <- model.matrix(~ ., df_assessment_X_nm)
  df_assessment_y <- df_assessment$medv
  
  #i番目のモデルに対して、検証セットにおける予測精度を確認します。
  model <- models_cv[[splits_cv$id[[i]]]]
  
  #λは事前に設定した候補(lambdas)分だけ検証します。
  for(lambda in lambdas){
    preds <- predict(
      object = model,
      Xnew = df_assessment_1X_nm,
      lambda = lambda
    )
    preds <- c(preds$fit)
    .rmse <- yardstick::rmse_vec(df_assessment_y, preds)
    
    #評価結果をデータフレームに記録します。
    results_cv <- dplyr::bind_rows(
      results_cv,
      tibble::tibble(id = splits_cv$id[[i]],
                     lambda = lambda,
                     rmse = .rmse))
  }
}

#各λでRMSEの平均値を計算し、昇順にソートします。
results_tmp <- results_cv |>
  dplyr::group_by(lambda) |>
  dplyr::summarise(mean_rmse = mean(rmse)) |>
  dplyr::arrange(mean_rmse)

#結果を表示します。
results_tmp |> head()

このλにおける予測精度は次のとおりです。

lambda_best <- results_tmp$lambda[[1]]

preds <- predict(
  object = model_genlasso,
  Xnew = model.matrix(~ ., df_test_X_nm),
  lambda = lambda_best
)
preds <- c(preds$fit)
.rmse <- round(yardstick::rmse_vec(df_test_y, preds), 3)
ggplot(mapping = aes(x = df_test_y, y = preds)) +
  geom_point() +
  ggtitle(paste("Lasso @", round(lambda_best, 3), " | RMSE:", .rmse))

1次元フューズドラッソ

ビニングを行う場合

フューズドラッソとは罰則項を\lambda \sum_{i=1}^{k-1} | \beta_{i+1} - \beta_{i} | としたもののことで、 隣接した説明変数同士の係数が同じ値に近づくような作用を与えます。

今回は例として、説明変数lstatと目的変数の関係に着目することとします。

横軸にlstat²、縦軸に目的変数をとって散布図を描くと次のとおりです。

lambda <- lambda_best
ggplot(mapping = aes(x = df_train_X_nm$lstat, y = df_train_y)) +
  geom_point() +
  ggtitle(paste("Lasso @", round(lambda, 3), " | RMSE:", .rmse))

この2変数の関係を、フューズドラッソにより要約して捉えることを考えましょう。

まず、lstatの値を離散化します。

値域をいくつかの領域（ビン）に区分し、 それぞれの領域に属しているときに1、そうでないときに0というダミー変数をビンの数だけ作成します。 なお、このような処理をビニングといいます。

X_raw <- df_train_X_nm$lstat
#領域を区切る点を指定します。
breaks <- seq(min(X_raw), max(X_raw), length.out = 21) #20等分

#各行に対する、ダミー変数化の処理をを定義します。
breaks_r <- c(breaks, +Inf)
fn_discretize <- function(x) {
  as.integer(x >= breaks_r[-length(breaks_r)] & x < breaks_r[-1])
}
#この処理を、全ての行で同時に行います。
X <- sapply(X_raw, fn_discretize)
#行と列が入れ替わってしまうため、転置します。
X <- t(X)

#もともとのlstatと、各ダミー変数の値を横に並べて表示します。
head(cbind(X_raw, X))

           X_raw                                          
[1,]  2.53339277 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
[2,]  0.09071772 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[3,] -0.76778399 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[4,] -0.07986406 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[5,] -0.72024481 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[6,] -0.88523374 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

同じことを後で評価用データにも行う必要があるので、関数化しておきます。

fn_getdisclstat <- function(df_X_nm){
  X <- df_X_nm$lstat
  X <- sapply(X, fn_discretize)
  t(X)
}

計画行列はこれで完成です。 フューズドラッソに対応する D はgetD1d関数で取得できるため、 次のように指定することで実行可能ですが……

model_fusedlasso1d <- genlasso(y = df_train_y, X = X, D = getD1d(ncol(X)))

フューズドラッソの場合はより実行が簡単で、かつ効率的なアルゴリズムを用いることが出来る fusedlasso1d関数が用意されています。

model_fusedlasso1d <- fusedlasso1d(y = df_train_y, X = X)

plot関数により、フューズドラッソの係数（要約した結果）をプロットすることができます。

plot(model_fusedlasso1d, lambda = lambda)

予測の方法などはgenlasso関数を用いた場合と同じです。

lambda_fused <- 0.027 * nrow(df_train)

preds <- predict(
  object = model_fusedlasso1d,
  Xnew = fn_getdisclstat(df_train_X_nm),
  lambda = lambda_fused
)
preds <- data.frame(lstat = df_train_X_nm$lstat, y = df_train_y, source = "actual") |>
  rbind(data.frame(lstat = df_train_X_nm$lstat, y = c(preds$fit), source = "model"))
ggplot(preds, mapping = aes(x = lstat, y = y, color = source)) +
  geom_point() +
  ggtitle(paste("1 Var Fused Lasso @", round(lambda_fused, 3)))

この1変数だけの予測モデルの評価用データにおける予測精度を確認すると次のとおりです。

preds <- predict(
  object = model_fusedlasso1d,
  Xnew = fn_getdisclstat(df_test_X_nm),
  lambda = lambda
)
preds <- c(preds$fit)
.rmse <- round(yardstick::rmse_vec(df_test_y, preds), 3)
ggplot(mapping = aes(x = df_test_y, y = preds)) +
  geom_point() +
  ggtitle(paste("1 Var Fused Lasso @", lambda, " | RMSE:", .rmse))

ビニングを行わない場合

ラッソ回帰の場合はモデルの係数の個数が非常に多くとも問題ないという特徴があるため、 ビニングしないでそのままのデータを用いることも考えられます。

Xを省略した場合はyの長さと同じ数だけの説明変数が用意されます。 このとき、「隣り合う説明変数」を正しく認識させるため、 yを元の説明変数の順序で並び替えておきます。

#見た目がそれなりに滑らかになるλに設定します。
lambda_fused2 <- 0.3 * nrow(df_train)
#説明変数の順序で並べ替えます。
y_sorted <- sort_by(df_train_y, X_raw)
#モデルを構築します。
model_fusedlasso1d_a <- fusedlasso1d(y = y_sorted)
#構築したモデルと元のデータを重ねてプロットします。
plot(model_fusedlasso1d_a, lambda = lambda_fused2)

引数posで位置情報を与えることも可能ですが、 フューズドラッソの場合は特に計算結果に影響しません。

#昇順という制限があるため、ソートしておきます。
pos <- sort(X_raw)
#重複なしという制限があるため、重複している箇所はわずかに数値をずらします。
pos <- pos + seq(0, 1e-7, length.out = length(X_raw))

model_fusedlasso1d_ap <- fusedlasso1d(y = y_sorted, pos = pos)
plot(model_fusedlasso1d_ap, lambda = lambda_fused2)

トレンドフィルター

罰則項を \lambda \sum_{i=1}^{k-2} | (\beta_{i+2} - \beta_{i+1}) - (\beta_{i+1} - \beta_{i}) | とすることで、 隣接した「説明変数の差」が同じ値に近づくような作用を与えることができます。 一般化ラッソ回帰問題にて D を次のように設定することでこの問題に帰着します。

D = \begin{bmatrix} 1 & -2 & 1 & 0 & \cdots & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & \cdots & 0 & 0 & 0 \\ & & \vdots & & \ddots & & \vdots & \\ 0 & 0 & 0 & 0 & \cdots & 1 & -2 & 1 \\ \end{bmatrix}

「なるべく傾きの変わらない（少ない数の）直線でつなぐ」ことで、長期的なトレンドを抽出することができます。

genlasso関数にて引数DをgetDtfPosSparse関数で与えることでも実行できますが、 より実行が簡単で、かつ効率的なアルゴリズムを用いることが出来る trendfilter関数が用意されています。

model_trendfilter1d<- trendfilter(y = y_sorted, ord = 1)
plot(model_trendfilter1d, lambda = lambda_fused2)

さらに高階の差分を考えることもできます。

つまり、罰則項を \lambda \sum_{i=1}^{k-3} | \beta_{i+3} - 3\beta_{i+2} + 3\beta_{i+1} - \beta_{i} | とすることで「なるべく少ない数の放物線でつなぐ」（2次のトレンドフィルター）、 罰則項を \lambda \sum_{i=1}^{k-4} | \beta_{i+4} - 4\beta_{i+3} + 6\beta_{i+2} - 4\beta_{i+1} + \beta_{i}| とすることで「なるべく少ない数の3次曲線でつなぐ」（3次のトレンドフィルター）という回帰問題とすることができます。

引数ordを変更することでこのような計算を行うことができます。 なお、ゼロの場合はフューズドラッソになります。

model_trendfilter2d<- trendfilter(y = y_sorted, ord = 2)
plot(model_trendfilter2d, lambda = lambda_fused2)

model_trendfilter3d<- trendfilter(y = y_sorted, ord = 3)
plot(model_trendfilter3d, lambda = lambda_fused2)

各点の位置posを与えることも可能です。 ordが1以上の場合は、罰則項の係数が隣り合う点の距離に応じて調整される（「1/距離」倍だけ補整される）ため、 結果が前述したものとは異なってきます。

ただし、今回のデータのように点が密集している箇所がある場合、あまり意味のある結果にはなりません。

model_trendfilter1d_p<- trendfilter(y = y_sorted, ord = 1, pos = pos)
plot(model_trendfilter1d_p, lambda = lambda_fused2)

model_trendfilter2d_p<- trendfilter(y = y_sorted, ord = 2, pos = pos)
plot(model_trendfilter2d_p, lambda = lambda_fused2)

model_trendfilter3d_p<- trendfilter(y = y_sorted, ord = 3, pos = pos)
plot(model_trendfilter3d_p, lambda = lambda_fused2)

なお、このtrendfilterに対しては 交差検証法（CV）でλをチューニングするcv.trendfilterという関数が存在しています。

model_trendfilter0d<- trendfilter(y = y_sorted, ord = 0)
cv <- cv.trendfilter(model_trendfilter0d)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

plot(model_trendfilter0d, lambda=cv$lambda.min, 
     main=paste0("lambda = ", cv$lambda.min))

2次元フューズドラッソ

2次元の数値データに対しては、 水平方向と垂直方向の両方の隣接関係に対する係数の差分を罰則とするよう D を設定することで、 例えば画像のノイズ除去のような処理を行うこともできます。

このような設定でのフューズドラッソ回帰を行う関数fusedlasso2dが用意されています。 実際に、この2次元フューズドラッソにより画像のノイズ除去を行ってみましょう。

まず、サンプル画像を読み込みます。

img_raw <- png::readPNG(system.file("img", "Rlogo.png", package="png"))
str(img_raw)

 num [1:76, 1:100, 1:4] 0 0 0 0 0 0 0 0 0 0 ...

サンプルデータは、縦76ピクセル、横100ピクセルの画像ファイルです。 各ピクセルはRGB（赤、緑、青）の3色の強さと、 不透明度を表すアルファチャンネルを合わせた4つのデータの組で表されるため、 縦ピクセル数×横ピクセル数×4 の3次元配列となってデータが格納されています。

これをRで描画するためにはrasterImage関数を用います。

fn_plotimg <- function(img, ...){
  img_width <- dim(img)[2]
  img_height <- dim(img)[1]
  
  par_old <- par(mar=c(1, 1, 1, 1)) #余白の調整
  if (img_width > img_height){
    plot(0:img_width, type='n', xlab = "", ylab = "", axes = FALSE, ...)
    rasterImage(img, 0, (img_width - img_height)/2, img_width, 
                img_height + (img_width - img_height)/2)
  }else{
    plot(0:img_height, type='n', xlab = "", ylab = "", axes = FALSE, ...)
    rasterImage(img, (img_height - img_width)/2, 0, 
                img_width + (img_height - img_width)/2, img_height)
  }
  par(par_old) #余白をもとに戻す
}
fn_plotimg(img_raw, main = "Raw Image")

これに、正規分布に基づくノイズを付加してみます。

set.seed(42)
img_noisy <- img_raw
img_noisy <- img_noisy + array(rnorm(length(img_noisy), sd = 0.3), dim(img_noisy))
img_noisy[img_noisy > 1] <- 1
img_noisy[img_noisy < 0] <- 0
fn_plotimg(img_noisy, main = "Noisy Image")

RGB、アルファチャンネルの4チャンネルそれぞれに対して、 2次元フューズドラッソによるノイズ除去を試みます。

ただし、回帰モデルの説明変数の数が「縦ピクセル数×横ピクセル数」となることから、 この画像全体に対して1つのモデルを構築しようとするとモデル構築に時間がかかってしまうため、 画像の一部を切り出してその部分にのみ処理を行うこととします。

λによってノイズと判断して除去される凸凹が変化することが分かります。 今回の例の場合、0.1～0.3あたりがちょうど良いようです。

#後で全体のノイズ除去を試みるにあたり、
#区画を size_window × size_window ピクセルずつに分けることとします。
size_window <- 15
#ただし、区画の境界で不自然な線が現れないよう、
#区画の右側と下側は size_pad ピクセル分広くとった領域でモデルを構築します。
size_pad <- 3

#モデルを格納するリストを用意します。
models_fusedlasso2d <- list()

c <- 1 #色のチャンネル
i <- 2 #Y軸方向のインデックス
j <- 3 #X軸方向のインデックス

#ノイズ除去を行う領域(上下左右の端)を設定します。
x1 <- (j - 1) * size_window + 1
y1 <- (i - 1) * size_window + 1
x2 <- min(x1 + size_window + size_pad - 1, dim(img_raw)[2])
y2 <- min(y1 + size_window + size_pad - 1, dim(img_raw)[1])

#cat("i = ",i, ", j = ",j, ", c = ",c,"\n")
#全体の画像から対象領域を取り出します。
img_tmp <- img_noisy[y1:y2,x1:x2,c]

models_fusedlasso2d[[i]] <- list()
models_fusedlasso2d[[i]][[j]] <- list()
#2次元フューズドラッソ回帰モデルを構築します。
models_fusedlasso2d[[i]][[j]][[c]] <- fusedlasso2d(y = img_tmp)

#2x2で画像を表示します。
par_old <- par(mfrow = c(2, 2))
fn_plotimg(img_raw[y1:y2,x1:x2,c], main = "Raw Image")
fn_plotimg(img_noisy[y1:y2,x1:x2,c], main = "Noisy Image")

#各λの結果を並べて表示します。
for(lambda in c(0.01, 0.03, 0.1, 0.3, 1, 3)){
  co <- coef(models_fusedlasso2d[[i]][[j]][[c]], lambda = lambda)
  img_processed_tmp <- array(co$beta, c(dim(img_tmp)[1], dim(img_tmp)[2]))
  fn_plotimg(img_processed_tmp, main = paste0("Processed Image, lambda = ",lambda))
}
#描画設定をもとに戻します。
par(par_old)

全体にこれを適用すると次のようになります。

今回の例のようなそれほど大きくない画像でも数十分程度の時間がかかるため、実行時は注意してください。

models_fusedlasso2d <- list()
ni <- ((dim(img_raw)[1]-size_pad-1) %/% size_window + 1)
nj <- ((dim(img_raw)[2]-size_pad-1) %/% size_window + 1)
for(i in 1:ni){
  models_fusedlasso2d[[i]] <- list()
  for(j in 1:nj){
    models_fusedlasso2d[[i]][[j]] <- list()
    
    x1 <- (j - 1) * size_window + 1
    y1 <- (i - 1) * size_window + 1
    x2 <- min(x1 + size_window + size_pad - 1, dim(img_raw)[2])
    y2 <- min(y1 + size_window + size_pad - 1, dim(img_raw)[1])
    for(c in 1:dim(img_raw)[3]){
      #cat("i = ",i, ", j = ",j, ", c = ",c,"\n")
      img_tmp <- img_noisy[y1:y2,x1:x2,c]
      models_fusedlasso2d[[i]][[j]][[c]] <- fusedlasso2d(y = img_tmp)
    }
  }
}

par_old <- par(mfrow = c(2, 2))
fn_plotimg(img_raw, main = "Raw Image")
fn_plotimg(img_noisy, main = "Noisy Image")
imgs_processed <- list()
for(lambda in c(0.01, 0.03, 0.1, 0.3, 1, 3)){
  img_processed <- array(1, dim = dim(img_raw))
  for(c in 1:dim(img_raw)[3]){#
    img_processed_tmp <- array(NA, dim = c(dim(img_raw)[1:2], ni*nj))
    for(i in 1:ni){
      for(j in 1:nj){
        x1 <- (j - 1) * size_window + 1
        y1 <- (i - 1) * size_window + 1
        x2 <- min(x1 + size_window + size_pad - 1, dim(img_raw)[2])
        y2 <- min(y1 + size_window + size_pad - 1, dim(img_raw)[1])
        co <- coef(models_fusedlasso2d[[i]][[j]][[c]], lambda = lambda)
        img_processed_tmp[y1:y2,x1:x2,nj*(i-1)+j] <- array(co$beta, c(y2-y1+1, x2-x1+1))
      }
    }
    #size_padによって広げた箇所は複数モデルの結果があるため、これらは単純平均します。
    img_tmp <- apply(img_processed_tmp, 1:2, mean, na.rm = T)
    img_tmp[is.na(img_tmp)] <- 1
    img_processed[,,c] <- img_tmp
  }
  imgs_processed[[as.character(lambda)]] <- img_processed
  fn_plotimg(imgs_processed[[as.character(lambda)]], main = paste0("Processed Image, lambda = ",lambda))
}
par(par_old)

他パッケージとの比較

glmnet

glmnet パッケージは正則化GLMのモデリングに用いられるパッケージです。

genlasso パッケージと比較すると次のような特徴があります。

• glmnet パッケージはGLMを取り扱うもので、よく使われるリンク関数や残差分布（ガウス分布以外）に対応しています。genlasso パッケージは通常の線形回帰問題のみです。
• genlasso パッケージ利用時は、標準化や切片項の追加を自分で行う必要がありましたが、glmnet パッケージでは自動で行えます。
• glmnet パッケージは各観測に対する重み付けにも対応しています。
• glmnet パッケージも genlasso パッケージ同様、一度の学習で全てのλに対する解を得ることができます。
• glmnet パッケージは、交差検証法（CV）によるλのチューニングにも対応しています。genlasso パッケージの場合、同機能があるのはトレンドフィルターのみです。
• glmnet パッケージは、ラッソ正則化のみならず、リッジ（2 - ノルムによる罰則項）やエラスティックネット（ラッソとリッジ両方の罰則項を持つもの）にも対応しています。一方、D\beta の形の罰則項には対応していません。

このように、予測モデリングに便利な機能が多数備わっているため、 glmnetパッケージが使用できる状況（罰則項が通常のラッソ回帰である場合）ではこちらの方が有用なケースが多いでしょう。

実際、 通常のラッソ回帰問題 節で述べた例は（λのチューニングも含めて） 次のとおり簡単に記述することができます。

library(glmnet)
#CVで最適なλを求めながらモデルを構築します。
set.seed(42)
model_cv.lasso <- glmnet::cv.glmnet(
  #matrix型である必要があるため変換します。
  x = as.matrix(df_train_X_nm),
  y = df_train_y,
  #正則化項のαを指定（1：ラッソ、0：リッジ）します。
  alpha = 1, 
  #CVの分割数を指定します。
  nfolds = 5 
)
#最適なλを抽出します。
lambda_glm <- model_cv.lasso$lambda.min |> round(3)

#予測精度の確認のため、評価用データでの予測値を取得します。
preds <- predict(
  object = model_cv.lasso,
  #matrix型である必要があるため変換します。
  newx = as.matrix(df_test_X_nm), 
  #λは引数sに与えます（lambdaではありません）。
  s = lambda_glm
) |> c()
.rmse <- round(yardstick::rmse_vec(df_test_y, preds), 3)
ggplot(mapping = aes(x = df_test_y, y = preds)) +
  geom_point() +
  ggtitle(paste("Lasso @", lambda_glm, " | RMSE:", .rmse))

#λごとのCVの結果をプロットします。
plot(model_cv.lasso)

#λ別の係数プロット（解パス）をプロットします。
plot(model_cv.lasso$glmnet.fit, "lambda")

#指定したλにおける係数を出力します。
#genlassoの係数とほぼ同等となっていることがわかります。
coef(model_cv.lasso$glmnet.fit, s = c(0.01, 0.1, 1, 10))

14 x 4 sparse Matrix of class "dgCMatrix"
                     s1         s2          s3       s4
(Intercept) 22.44828496 22.4482850 22.44828496 22.44828
crim        -1.07290435 -0.7884986  .           .      
zn           0.81966961  0.4262189  .           .      
indus       -0.03945084 -0.1777379  .           .      
chas         0.47170645  0.4563839  .           .      
nox         -2.10738128 -1.6769990  .           .      
rm           2.50046418  2.6201820  2.54685806  .      
age         -0.16251918  .          .           .      
dis         -3.23779431 -2.4650921  .           .      
rad          2.61941760  1.2670168  .           .      
tax         -1.72933402 -0.5630096  .           .      
ptratio     -2.15099633 -2.0243315 -1.22146843  .      
black        0.68711370  0.6314864  0.05446536  .      
lstat       -3.84813288 -3.8995978 -3.66282065  .      

なお、glmnet パッケージで最小化される関数は次のようなものです。

\frac1n\sum_{i=1}^nw_il \left(y_i, \beta_0 + \sum_{j=1}^{k} \left( x_{ij} \cdot \beta_j \right) \right) + \lambda \left[ \frac{(1-\alpha)}2 \| \beta \| _2^2 + \alpha \| \beta \|_1 \right]

ただし、\beta_0 はモデルの切片項（罰則項のノルムの計算には含まれません）、 l(y_i, \eta_i) は i 番目の観測に対する負の対数尤度関数、 w = (w_i)_{i=1}^n は観測の重みを表すベクトル、 \alpha \in [0, 1] は罰則項の形状を表すハイパーパラメータです。

リンク関数として恒等関数、残差分布としてガウス分布を選ぶと、 l(y_i, \eta_i) = \frac12 (y_i - \eta_i)^2 となり、 また切片項 \beta_0 = 0 、重み w を全て1、\alpha = 1 （ラッソ回帰）とすることで genlasso が最小化する関数と同等の形式になります。

\frac1{2n}\| y - X\beta \|_2^2 + \lambda\| \beta \|_1

ただし、第1項に係数 1/n がかかっている点のみ genlasso のものとは異なっています。 これは、第1項が表す残差が観測数に比例して増加する点を補整するためです。

このため、解として得られる λ のスケールが異なることに注意してください。 具体的には、genlasso 側の λ が観測数 n に比例して大きくなります ³。

aglm

aglm とは正則化 GLM の変種であるAGLM（Accurate Generalized Linear Model） Suguru Fujita et al. (2020) のモデリングを実装したパッケージです。

AGLM は GLM とよく似たモデルですが、 連続変数に対しては ビニングを行う場合 で述べたような ビニングを伴うフューズドラッソと同等の変換を行います。 GLM のような加法的モデルであるという性質を維持しながら、 説明変数と目的変数の間の非線形な関係を表現することが可能になっています。

ここでダミー変数化の際に順序付きダミー変数（閾値以上なら 1 となるダミー変数）を利用することで、 フューズドラッソが通常の（罰則項を \lambda \| \beta \| とする）ラッソ回帰に帰着できるという性質を用いており、 ダミー変数化の処理を行った後の内部の計算処理に glmnet を用いているのが特徴です。 そのため、前述した glmnet の特徴（利点）がそのまま継承されています。

次のようにすることで、ビニングを行う場合 の例を再現するモデルを構築することができます。

library(aglm)
model_aglm_1v <- aglm(
  x = df_train_X_nm$lstat,
  #ビンの端点をfusedlasso1dで使用したものと同じにします。
  bins_list = list(X1 = breaks), 
  y = df_train_y,
  add_linear_columns = FALSE,
  #glmnetで説明変数を標準化しないようにします。
  standardize = F, 
  #OD化する際に離散化する（線形補間しない）ようにします。
  OD_type_of_quantitatives = "J", 
  alpha = 1 # 正則化項のαパラメーターを指定（1：ラッソ、0：リッジ）
)
#λはgenlassoの例から流用します。
lambda_aglm <- lambda_best / nrow(df_train_X)

#予測精度の確認のため、評価用データでの予測値を取得します。
preds <- predict(
  object = model_aglm_1v,
  newx = df_test_X_nm$lstat,
  s = lambda_aglm
) |> c()
.rmse <- round(yardstick::rmse_vec(df_test_y, preds), 3)
ggplot(mapping = aes(x = df_test_y, y = preds)) +
  geom_point() +
  ggtitle(paste("1 Var Fused Lasso by AGLM @", round(lambda_aglm, 3),
                " | RMSE:", .rmse))

#関数の形をプロットします。
par(mar = c(4, 4, 1, 1)) # プロット領域中の余白を調整
plot(model_aglm_1v, verbose = FALSE, vars = NULL, s = lambda_aglm, ask = FALSE)

これは説明変数 lstat にのみ着目した1変数のモデルでしたが、 全変数を説明変数としてモデルを構築すると次のとおりです⁴。

#CVで最適なλを求めながらモデル構築
set.seed(42)
model_cv.aglm <- cv.aglm(
  x = df_train_X,
  y = df_train_y,
  add_linear_columns = FALSE,
  alpha = 1, # 正則化項のαパラメーターを指定（1：ラッソ、0：リッジ）
  nfolds = 5 # クロスバリデーションの分割数を指定
)
lambda_aglm_best <- model_cv.aglm@lambda.min |> round(3)

非線形な関係を捉えることができているため、 GLMと比べて予測精度が大幅に向上していることがわかります。

#予測精度の確認
preds <- predict(
  object = model_cv.aglm,
  newx = df_test_X,
  s = lambda_aglm_best
) |> c() # 行列形式の出力をベクトルに変換
.rmse <- round(yardstick::rmse_vec(df_test_y, preds), 3)
ggplot(mapping = aes(x = df_test_y, y = preds)) +
  geom_point() +
  ggtitle(paste("AGLM@", lambda_aglm_best, " | RMSE:", .rmse))

予測関数は説明変数の1変数関数の和で表現されているので、 説明変数ごとに予想関数を分けてプロットすることで 予測関数の形を把握することも可能です。

par(mar = c(4, 4, 1, 1)) # プロット領域中の余白を調整
plot(model_cv.aglm, verbose = FALSE, vars = NULL, 
     s = lambda_aglm_best, ask = FALSE)

このように AGLM は予測精度と解釈可能性をある程度両立したモデルであり、 さらに glmnet パッケージの利点を多く継承していることから、 genlasso パッケージよりもこちらを利用したほうがよいケースも多いでしょう。

参考文献

Ryan J. Tibshirani, and Jonathan Taylor. 2011. “The Solution Path of the Generalized Lasso.” https://www.stat.cmu.edu/~ryantibs/papers/genlasso.pdf.

Suguru Fujita, Toyoto Tanaka, Kenji Kondo, and Hirokazu Iwasawa. 2020. “A Hybrid Modeling Method of GLM and Data Science Techniques.” https://www.institutdesactuaires.com/global/gene/link.php?doc_id=16273&fg=1.

Taylor B. Arnold, and Ryan J. Tibshirani. 2016. “Efficient Implementations of the Generalized Lasso Dual Path Algorithm.” https://www.stat.cmu.edu/~ryantibs/papers/fastgl.pdf.

Taylor B. Arnold, and Ryan Tibshirani. 2019. “Introduction to the Genlasso Package.” https://ryantibs.r-universe.dev/genlasso/doc/article.pdf.

岩沢 宏和, and 平松 雄司. 2019. 入門 Rによる予測モデリング. 東京図書.

藤田 卓, 田中 豊人, and 岩沢 宏和. 2019. “アクチュアリー実務のためのデータサイエンスの技術を用いたglmの拡張.” https://www.jarip.org/publication/risk_and_insurance/pdf/RI_v15_045.pdf.

---

1. ここでは切片項を省略して表示しています。

2. 以下断りが無ければ、通常のラッソ回帰問題 で標準化したものです。

3. 本稿ではこれを踏まえて、定数×観測数 の形で λ を指定している箇所が存在します。

4. 直上の例での標準化やビンの設定などは genlasso パッケージの結果を再現するために設定したものです。 以下の例では再現の必要が無いため、aglm パッケージのデフォルト値に戻しています。