scikit-learnで求めたLassoの解パスとCVの結果を可視化する

nayu.t.s

2022.09.25

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

データアナリティクス事業本部の鈴木です。

Lassoは、例えばscikit-learnを使うと簡単にクロスバリデーション（以降CV）で正則化パラメータとそのときの回帰係数を推定できます。一方で、単にCVで最適なモデルを求めるだけではなく、解パス図を描いて回帰係数の変化を詳しくみたり、CVで求めた値を使ってより簡単なモデルを選択したりしたいと思ったので、scikit-learnの関数を使って試してみました。

Lassoとは

Lassoは線形回帰モデルの回帰係数の推定法の一つで、係数の推定時にいくつかの係数を0にすることで変数選択を行える手法です。

以下のように回帰係数のL1ノルムを正則化項として誤差二乗和に付けた関数を、回帰係数について最小化することで回帰係数を推定します。

lasso

ただし、L1ノルムは以下です。

L1ノルム

例えば2変数の場合、以下の斜線部に誤差二乗和部分（第1項目）の解β^LSが入る場合にそのβの対応する成分は0になります。※参考1のP.17

回帰係数のいずれかが0になる領域

中央のひし形は正則化付きの解βの実行可能領域に対応します。ひし形は正則化パラメータαによって大小し、この内部にβが値を取ります。斜線部に誤差二乗和部分の解β^LSが入る場合には、ひし形の中で取れる解βが最もその解に近づけるのがひし形の頂点なので、その頂点に対応する成分は0になるという解釈ができます。

誤差二乗和の項も正則化項も共に凸なので、全体も凸関数となり、最適化すれば解βを求めることができます。

最適化については、Lassoの解を求めるための方法がいろいろありますが、scikit-learnのAPIを使うことでそこはお任せすることができます。

実行環境・ノートブック

今回はGoogle Colaboratory環境で実行しました。

ライブラリのバージョンは以下です。

numpy: 1.21.6
scipy: 1.7.3
scikit-learn: 1.0.2

ノートブックは以下のGithubレポジトリに格納してあります。

やってみる

上記レポジトリのJupyter Notebookに全体の流れを記載してあるので、ここではポイントのみを記載します。

解パス図の描写

lasso_path関数を使って、αとそのときの回帰係数を得ることができます。係数が得られれば、係数ごとに見やすいよう色を変えるなどして、αごとにプロットします。

このように横軸に正則化パラメータ、縦軸に各回帰係数を取ってプロットすることで、αと回帰係数の関係を表した図を解パス図と言います。

以下のコードで描写してみます。

# 以下のコードを修正しました。
# https://scikit-learn.org/stable/auto_examples/linear_model/plot_lasso_coordinate_descent_path.html#sphx-glr-auto-examples-linear-model-plot-lasso-coordinate-descent-path-py
# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
# License: BSD 3 clause

# データをlasso_pathに渡して解パスを求める
alphas_lasso, coefs_lasso, _ = lasso_path(X, y, eps=eps)

# 解パスを描写する
plt.figure(1)
colors = cycle(["b", "r", "g", "c", "k"])
neg_log_alphas_lasso = -np.log10(alphas_lasso)
for coef_l, c in zip(coefs_lasso, colors):
    l1 = plt.plot(neg_log_alphas_lasso, coef_l, c=c)

plt.xlabel("-Log(alpha)")
plt.ylabel("coefficients")
plt.title("Lasso Paths")
l1[-1].set_label("Lasso")
plt.legend()
plt.axis("tight")

今回のデータだと、以下のようになりました。今は横軸が正則化パラメータαの対数値に-1を掛けたものなので、左に行くほど、つまりαが大きくなるほど0になる回帰係数が増えていく様子が分かります。（横軸は結果が比較しやすいよう、参考にしたドキュメントに合わせただけです。）

解パスの描写例

今回は色分けだけにしていますが、回帰係数ごとに対応する変数名が分かるようにすると、αによってどの回帰係数がどのように遷移するか分かり、そのデータに対する変数選択がLassoによってどのように行われるかよりよく分かります。

CVの結果の確認

Lassoの正則化パラメータαの選び方の一つとして、CVの結果で最もスコアが小さいものを採用する方法が考えられます。scikit-learnの場合、LassoCVが提供されています。

LassoCVのfit関数で最適化することにより、最も良いモデルを選択できるのでやってみます。

# 以下のコードを修正しました。
# https://scikit-learn.org/stable/auto_examples/linear_model/plot_lasso_model_selection.html#sphx-glr-auto-examples-linear-model-plot-lasso-model-selection-py
# Author: Olivier Grisel
#         Gael Varoquaux
#         Alexandre Gramfort
#         Guillaume Lemaitre
# License: BSD 3 clause

# データの追加
rng = np.random.RandomState(42)
n_random_features = 14
X_random = pd.DataFrame(
    rng.randn(X.shape[0], n_random_features),
    columns=[f"random_{i:02d}" for i in range(n_random_features)],
)
X = pd.concat([X, X_random], axis=1)

# CVの実行
lasso = LassoCV(cv=20, random_state=0).fit(X, y)

# 結果のグラフ化
plt.semilogx(lasso.alphas_, lasso.mse_path_, linestyle=":")
plt.plot(
    lasso.alphas_,
    lasso.mse_path_.mean(axis=-1),
    color="black",
    label="Average across the folds",
    linewidth=2,
)

plt.axvline(lasso.alpha_, linestyle="--", color="black", label="Best alpha")

plt.xlabel("alpha")
plt.ylabel("MSE")
plt.title("MSE for each alpha")

plt.legend()

図で確認すると、縦の破線に相当するαが選べました。

もう少し選択の根拠を調べてみると、MSEの平均値が最小になるαを使っていることが分かります。

print(f"Best Alpha (from model): {lasso.alpha_}")
# Best Alpha (from model): 0.6856153997937577
print(f"Best Alpha (from MSE path): {lasso.alphas_[np.argmin(lasso.mse_path_.mean(axis=-1),)]}")
# Best Alpha (from MSE path): 0.6856153997937577

ほかにも、MSEの標準誤差も加味して、1標準誤差ルール(※参考2のP.22)という方法で正則化パラメータを選択することもできます。

まず各αにおけるMSEの標準誤差をエラーバーで表現してみます。

# 標準誤差の計算
sem_err = [ sem(lasso.mse_path_[i,:]) for i in range(lasso.mse_path_.shape[0])]

# グラフの作成
plt.figure(figsize=(10,10))
plt.semilogx(lasso.alphas_, lasso.mse_path_, linestyle=":")
plt.plot(
    lasso.alphas_,
    lasso.mse_path_.mean(axis=-1),
    color="black",
    label="Average across the folds",
    linewidth=1,
)
plt.errorbar(lasso.alphas_, lasso.mse_path_.mean(axis=-1), yerr = sem_err, capsize=2, fmt='o', ecolor='k', ms=5, mfc='None', mec='k')

plt.axvline(lasso.alpha_, linestyle="--", color="green", label="alpha CV")

plt.ylim([1300, 7000])
plt.xlabel("alpha")
plt.ylabel("Standard Error of MSE")
plt.title("MSE for each alpha with standard error")

plt.legend()

エラーバー付きのαごとのMSE

図示する過程でMSEの標準誤差が計算できたので、これを使って1標準誤差ルールでの正則化パラメータαを選択してみます。

1標準誤差ルールは、最小のMSEの平均値の値とそのときの標準誤差を足した値を超えない範囲で、最大のMSEの平均値を取るαを選びます。

# CVで選んだ最適なalphaでのMSEの平均値とその標準誤差を調べる
best_alpha_idx = np.argmin(lasso.mse_path_.mean(axis=-1),)
se0 = sem_err[best_alpha_idx]
cv0 = lasso.mse_path_.mean(axis=-1)[best_alpha_idx]

# CV0+se0に収まる範囲の最大のMSE（の平均値）を取るalphaを求める
scores = lasso.mse_path_.mean(axis=-1)
scores = [score if score < cv0+se0 else -np.inf for score in scores]
alpha_one_standard_error_rule = lasso.alphas_[np.argmax(scores)]
print(f"Alpha: {alpha_one_standard_error_rule}")
# Alpha: 5.963133101500133

先ほどのプロットでαがどのあたりか図示すると、青色の縦線の位置になります。

1標準誤差ルールで求めたα

計算上当たり前ではありますが、より大きなαを選択できました。先ほど作った解パスと合わせて考えると、αが大きいためより多くの変数が0になり、変数の数が少ない、より簡単なモデルになっていることが想像できます。

CVを使った選択は、特に最近だと多くの機械学習手法で使われるため多くの人にも分かりやすく良い方法ですが、変数選択に関する一致性が保証されていないことが指摘されており、変数選択の手法の側面があるLassoを使う上では念頭に置いておくことが必要そうです。一致性を気にする場合は拡張BICなどの情報量基準を使うとよいですが、情報量基準を使った方法はscikit-learnではLassoLarsICとして実装されていました。

以下に例もあったため、興味がある方はご参照ください。