Lassoの理論と実装 -スパースな解の推定アルゴリズム-

はじめに

前回は、「重回帰モデルの理論と実装 -なぜ正則化が必要か-」ということで、L2正則化について最小二乗推定量を求めるところまでやりました。今回はその続きということで、L1正則化について取り上げたいと思います。
コードはこちらです。

正則化

重回帰モデル式(詳細は前回の記事を参照してください)、

y=Xβ+ϵ\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}

があったときに、LqL_q正則化を導入した二乗誤差関数は、

Sλ(β)=yXβ2+λβq=yXβ2+λi=0pβiq()\begin{aligned} S_{\lambda}(\boldsymbol{\beta}) & = ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda||\boldsymbol{\beta}||_q \\ & = ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda \sum_{i=0}^{p}{|\beta_i|^q} \cdots (*) \\ \end{aligned}

となります。

確認ですが、

  • y\boldsymbol{y}:n次元の観測ベクトル
  • XX:n×(p+1)次元の計画行列(design matrix)
  • β\boldsymbol{\beta}:(p+1)次元の回帰係数ベクトル (これを求めたい。)
  • ϵ\boldsymbol{\epsilon}:n次元の誤差ベクトル

となります。

L2正則化(Ridge)

特に、()(*)について、q=2q=2のものをL2正則化(Ridge)といいます。

Sλ(β)=yXβ2+λβ2S_{\lambda}(\boldsymbol{\beta}) = ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda||\boldsymbol{\beta}||^2

これに関しては、前回の記事に書いた通り、損失関数Sλ(β)S_{\lambda}(\boldsymbol{\beta})β\boldsymbol{\beta}について偏微分した値が00になるようなβ\boldsymbol{\beta}を計算すると、

βridge=arg minβRp[yXβ2+λβ2]=(XTX+λIp+1)1XTy\begin{aligned} \boldsymbol{\beta}_{\text{ridge}} & = \argmin_{\boldsymbol{\beta} \in \mathcal{R^p}} \biggl[ ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda ||\boldsymbol{\beta}||^2 \biggr] \\ & = (X^{\mathsf{T}}X + \lambda I_{p+1})^{-1}X^{\mathsf{T}}\boldsymbol{y} \end{aligned}

と求めることができました。

L1正則化(Lasso)

()(*)について、q=1q=1のものをL1正則化(Lasso)[Tibshirani, 1996]といいます。LassoはLeast absolute shrinkage and selection operatorの略らしいですね。

Sλ(β)=yXβ2+λβ1\begin{aligned} S_{\lambda}(\boldsymbol{\beta}) = ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda||\boldsymbol{\beta}||_1 \end{aligned}

L2正則化のときと同様に、β\boldsymbol{\beta}で偏微分して。。。推定量を求めたいところですが、L1の誤差関数の正則化項(i.e.,λβi.e., \lambda |\boldsymbol{\beta}|)がβ\boldsymbol{\beta}で偏微分不可能な点を含むため、L2正則化のときのようには推定量を求めることができません。そこで今回はCD(Coordinate Descent)[J Friedman et al., 2007;2010]というアルゴリズムを用いてβlasso\boldsymbol{\beta}_{\text{lasso}}を推定したいと思います。

βlasso=arg minβRp [yXβ2+λβ1]\begin{aligned} \boldsymbol{\beta}_{\text{lasso}} & = \argmin_{\boldsymbol{\beta} \in \mathcal{R^p}} \biggl[ ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda ||\boldsymbol{\beta}||_1 \biggr] \end{aligned}

なぜL1正則化か

その前に、L1正則化をすると何が嬉しいのか、どういった回帰係数(β\boldsymbol{\beta})が得られるのか確認してましょう。 L1正則化を行うと、推定したパラメータ(今回は回帰係数であるβ\boldsymbol{\beta})の一部が0になり、スパースな解が得られます。特徴量がものすごくたくさんあるとき、つまり重回帰モデルおける説明変数がものすごくたくさんあるとき、それらの一部を0にすることによって、変数選択(特徴選択)してくれていることになります。 変数選択する指標として、

  • AIC
  • BIC
  • 各変数ごとの有意性検定

などが挙げられます。従来ではモデルの推定と変数選択は別々に行っていました。 Lassoではモデルの推定と変数選択を同時に行ってくれます。

次にどうしてLassoはスパースな解が得られるかということですが、次のサイトが参考になるかと思われます。

Coordinate Descentによる推定

アルゴリズム

CD(Coordinate Descent)については、

が非常に参考になるかと思われます。また、他の推定アルゴリズムとして有名なものは、

などがあります。

CDでは、各パラメータ(i.e., β1,β2,...,βpi.e.,\ \beta_1, \beta_2,...,\beta_p)毎に誤差関数を微分して更新式を得て、それを用いて更新を繰り返し行うことにより、収束した最適な推定値得ることができます。

Sλ(β)=12nyXβ2+λβ1\begin{aligned} S_{\lambda}(\boldsymbol{\beta}) & =\frac{1}{2n} ||\boldsymbol{y}-X\boldsymbol{\beta}||^2 + \lambda||\boldsymbol{\beta}||_1 \end{aligned}

今回は誤差関数を上のような式として、βj\beta_jで微分して更新式を得ます。 βj0\beta_j \neq 0のとき、βj\beta_jに関して微分し、=0=0とおいて、

1n[X:,jTX:,jβj+X:,jT(X:,jβjy)]+λsign(βj)=0\begin{aligned} \frac{1}{n} \biggl[ X_{:,j}^{\mathsf{T}}X_{:,j}\beta_j +X_{:,j}^{\mathsf{T}}(X_{:,-j}\boldsymbol{\beta}_{-j} - \boldsymbol{y}) \biggr] + \lambda \text{sign}(\beta_j) = 0 \end{aligned}

βj\beta_jについて解くと、

βj=1X:,jTX:,jS(X:,jT(yX:,jβj),λn)()\begin{aligned} \beta_j = \frac{1}{X_{:,j}^{\mathsf{T}}X_{:,j}} S\biggl(X_{:,j}^{\mathsf{T}}(\boldsymbol{y} - X_{:,-j}\boldsymbol{\beta}_{-j}), \lambda n \biggr) \cdots (**) \end{aligned}

ここで、SSsoft-thresholding operatorと呼ばれ、

S(x,λ)={xλ(x>λ)0(xλ)x+λ(x<λ)\begin{aligned} S(x, \lambda) = \left\{ \begin{array}{l} x - \lambda & (x > \lambda) \\ 0 & (|x| \leq \lambda) \\ x + \lambda & (x < -\lambda) \end{array} \right. \end{aligned}

という式で表されます。具体的にはλ=1\lambda=1のとき、横軸xx、縦軸S(x,λ)S(x, \lambda)のグラフは次のようになります。

また、切片(intercept)であるβ0\beta_0については、通常は正則化を考慮しません。2乗誤差をβ0\beta_0について微分し、求める事ができ、

β0=i=1n(yiXi,:β)n()\begin{aligned} \beta_0 = \frac{\sum_{i=1}^{n}( y_i-X_{i,:}\boldsymbol{\beta})}{n}\cdots (***) \end{aligned}

となります。 あとはβi\beta_iの更新を1,,p1,\cdots, pについて繰り返し行えばいいだけです! 更新アルゴリズムの擬似コードを示します。

LassoにおけるCoordinate Descent
01 β\boldsymbol{\beta}の初期化、λ\lambdaの固定
02 if 切片がある:
03   ()(***)式でβ0\beta_0更新。
04 while 収束条件を満たさない:
05   for jj in 1,,p1,\cdots, p:
06     ()(**)式でβj\beta_jを更新。
07   if 切片がある:
08     ()(***)式でβ0\beta_0更新。

Pythonによる実装

2016/01/11追記:コードをGithubのsatopirka/Lassoに上げました。

以下の実装では、説明変数同士を同じ尺度で評価するためにXXの各説明変数毎に基準化を行っています。 基準化(normalization)はいうまでもなく、平均値(μ\mu)を引いて標準偏差(σ\sigma)で割ることにより、平均値0分散1にすることです。

zi=xiμσz_i = \frac{x_i - \mu}{\sigma}

ちょっとだけscikit-learn風に書いてみました。

import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin
from copy import deepcopy

class Lasso(BaseEstimator, RegressorMixin):
  def __init__(self, alpha=1.0, max_iter=1000, fit_intercept=True):
    self.alpha = alpha # 正則化項の係数
    self.max_iter = max_iter # 繰り返しの回数
    self.fit_intercept = fit_intercept # 切片(i.e., \beta_0)を用いるか
    self.coef_ = None # 回帰係数(i.e., \beta)保存用変数
    self.intercept_ = None # 切片保存用変数

  def _soft_thresholding_operator(self, x, lambda_):
    if x > 0 and lambda_ < abs(x):
      return x - lambda_
    elif x < 0 and lambda_ < abs(x):
      return x + lambda_
    else:
      return 0

  def fit(self, X, y):
    if self.fit_intercept:
      X = np.column_stack((np.ones(len(X)),X))

    beta = np.zeros(X.shape[1])
    if self.fit_intercept:
      beta[0] = np.sum(y - np.dot(X[:, 1:], beta[1:]))/(X.shape[0])
    
    for iteration in range(self.max_iter):
      start = 1 if self.fit_intercept else 0
      for j in range(start, len(beta)):
        tmp_beta = deepcopy(beta)
        tmp_beta[j] = 0.0
        r_j = y - np.dot(X, tmp_beta)
        arg1 = np.dot(X[:, j], r_j)
        arg2 = self.alpha*X.shape[0]

        beta[j] = self._soft_thresholding_operator(arg1, arg2)/(X[:, j]**2).sum()

        if self.fit_intercept:
          beta[0] = np.sum(y - np.dot(X[:, 1:], beta[1:]))/(X.shape[0])

    if self.fit_intercept:
      self.intercept_ = beta[0]
      self.coef_ = beta[1:]
    else:
      self.coef_ = beta

    return self

  def predict(self, X):
    y = np.dot(X, self.coef_)
    if self.fit_intercept:
      y += self.intercept_*np.ones(len(y))
    return y

実際に使ってみます。前回の記事と同様にThe Boston Housing Datasetを使用します。またscikit-learn実装されているLassoと比較してみました。

  • コード
import pandas as pd

df = pd.read_csv("Boston.csv", index_col=0)
y = df.iloc[:,  13].values
df = (df - df.mean())/df.std() # 基準化
X = df.iloc[:, :13].values

model = Lasso(alpha=1.0, max_iter=1000)
model.fit(X, y)
print("Lasso")
print(model.intercept_)
print(model.coef_)

print("")
import sklearn.linear_model as lm

model = lm.Lasso(alpha=1.0, max_iter=1000, tol=0.0) # tol=0.0で収束判定なし(上の実装とほぼ同条件?)
model.fit(X, y)
print("")
print("scikit-learnのLasso")
print(model.intercept_)
print(model.coef_)
  • 出力
Lasso
22.5328063241
[ 0.          0.          0.          0.          0.          2.71517992
  0.          0.          0.          0.         -1.34423287  0.18020715
 -3.54700664]

C:\Python27_64\lib\site-packages\sklearn\linear_model\coordinate_descent.py:466: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations
  ConvergenceWarning)

scikit-learnのLasso
22.5328063241
[-0.          0.         -0.          0.         -0.          2.71517992
 -0.         -0.         -0.         -0.         -1.34423287  0.18020715
 -3.54700664]

一部の回帰係数が0となり、スパースな解が得られました。また、今回実装した結果はscikit-learnの結果と同じで、ちゃんと実装出来ていることが確認できました。

最後に、正則化パラメータλ\lambdaの大きさを変えたときに、得られる係数の値について見ていきたいと思います。λ\lambdaを大きくすればするほど、正則化項の値が大きくなります。つまり、回帰係数β\boldsymbol{\beta}のうち0となるものが増えていきます。逆にλ\lambdaを小さくしていけば、通常の重回帰モデルと同じくなるので、回帰係数β\boldsymbol{\beta}のうち0となるものが減っていきます。λ\lambdaを変化させたときに得られる係数をプロットしてみたのが、以下の図です。

うっ色が被っててわかりにくい。。。

最後に

今回は重回帰モデルにおけるL1正則化、Lassoについて扱いました。

参考


Written by@satopirka
A software engineer.

GitHubTwitter

© 2021 Minato Sato. All Rights Reserved