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

October 22, 2017

はじめに

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

正則化

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

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

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

\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}

となります。

確認ですが、

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

となります。

L2正則化(Ridge)

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

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

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

\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=1$ のものをL1正則化(Lasso)[Tibshirani, 1996]といいます。LassoはLeast absolute shrinkage and selection operatorの略らしいですね。

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

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

\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)については、

A Coordinate Descent Algorithm for the Lasso Problem - Jocelyn T. Chi
lecture 24: coordinate descent - YouTube (カーネギーメロン大学の凸最適化の講義動画)

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

Least angle regression(LARS)[Efron et al, 2004]

などがあります。

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

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

今回は誤差関数を上のような式として、 $\beta_j$ で微分して更新式を得ます。 $\beta_j \neq 0$ のとき、 $\beta_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}

$\beta_j$ について解くと、

\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}

ここで、 $S$ はsoft-thresholding operatorと呼ばれ、

\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}

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

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

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

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

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

Pythonによる実装

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

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

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$ を変化させたときに得られる係数をプロットしてみたのが、以下の図です。