はじめに

重回帰分析の発展として正則化について勉強しました。今回はラッソ回帰(L1正則化)のスクラッチ実装に挑戦しています。前回ラッソの理論についてまとめた記事は下記になりますので、興味のある方はご確認いただければと思います。

ラッソ回帰(L1正則化)を理解する(理論編)

参考

ラッソ回帰(L1正則化)の理解に当たって下記を参考にさせていただきました。

ラッソ回帰(L1正則化)の理論の復習

ラッソ回帰(L1正則化)の概要

$$L = \dfrac {1}{2}(\boldsymbol{y} - X\boldsymbol{w})^T (\boldsymbol{y} - X\boldsymbol{w}) + \lambda|| \boldsymbol{w} ||_{1}$$

上記がラッソ回帰(L1正則化)の損失関数の式になります。重回帰分析の損失関数に正則化項 $\lambda|| \boldsymbol{w} ||_{1} $を付け足した形になっています。ラッソ回帰(L1正則化)では上記のように重み$\boldsymbol{w}$のL1ノルムを加えることで正則化行います。

L1ノルムとは何か

ベクトル成分の絶対値の和(マンハッタン距離と呼ばれる)がL1ノルムです。ノルムは「大きさ」を表す指標で他にL1ノルムやL∞ノルムなどが使われます。

L1正則化の効果

L1正則化項を損失関数に加えることで重み$\boldsymbol{w}$の値を0へと誘導します。
その結果、余計な説明変数(特徴量)をモデルから省く次元圧縮的な効果があります。

ラッソ回帰(L1正則化)で解くべきタスク

$$L = \dfrac {1}{2}(\boldsymbol{y} - X\boldsymbol{w})^T (\boldsymbol{y} - X\boldsymbol{w}) + \lambda|| \boldsymbol{w} ||_{1}$$

上記損失関数を最小化するのがラッソ回帰の解くべきタスクでした。ラッソ回帰は正則化項に絶対値を含んでいるため、場合分けをしながら右微分や左微分といった方法を駆使して最適な重みの導出を目指します。

最適な重みを導出するためのラッソ回帰のアルゴリズムは複数ありますが、今回実装するのは座標降下法（coordinate descent）という手法です。座標降下法は重みを1つずつ更新する操作を繰り返して最適解に近似させる方法です。複数の重み$w_{0},w_{1},w_{2},...w_{n}$があるとすると、一度に全ての重みを更新するのではなく重みを1つずつ順番に更新します。

また前回の記事では上記数式を解くべきタスクとしていましたが、今回は損失関数を下記のようにおきます。

$$L = \dfrac {1}{2n}(\boldsymbol{y} - X\boldsymbol{w})^T (\boldsymbol{y} - X\boldsymbol{w}) + \lambda|| \boldsymbol{w} ||_{1}$$

1/2を1/2nに変更しています。スクラッチ実装した結果とsklearの結果を照らし合わせて検証するため、sklearで実装されている方法に合わせています。1/2nに変更することによって、重みのパラメータの値をより$0$に近づけていく効果があるようです。

ラッソ回帰(L1正則化)の重みの更新式

ラッソ回帰の各パラメータの更新式はこのようになっていました。
この式をもとにラッソ回帰の実装を進めていきます。
詳細については前回の記事ご参照ください。

バイアスの更新式

w_{0}=\frac{1}{j}\sum_{k=1}^{j}(y_{k}-\sum_{l=1}^{i}w_{k}x_{kl})

残りの重みのパラメータの更新式

右偏微分した際の更新式

\begin{eqnarray}

w_n^+ &=& \frac{\sum_{k=1}^{j}x_{kn}(y_{k} - w_{0} - \sum_{l \neq    n}^{i}w_{l}x_{kl})- n\lambda}{\sum_{k=1}^{j}x_{kn}^{2}}

\end{eqnarray}

左偏微分した際の更新式

\begin{eqnarray}

w_n^- &=& \frac{\sum_{k=1}^{j}x_{kn}(y_{k} - w_{0} - \sum_{l \neq    n}^{i}w_{l}x_{kl})+ n\lambda}{\sum_{k=1}^{j}x_{kn}^{2}}

\end{eqnarray}

ラッソ回帰の重み更新時の制約条件

下記のような時は$w_{n}$は初期値$0$のまま留まる。

\begin{eqnarray}

-n\lambda ≦ \sum_{k=1}^{j}x_{kn}(y_{k} - w_{0} - \sum_{l \neq  1}^{i}w_{l}x_{kl}) ≦ n\lambda

\end{eqnarray}

ラッソ回帰のスクラッチ実装

ラッソ回帰の実装

下記がラッソ回帰をスクラッチ実装してみたものです。
主に下記本と記事を参考にさせていただきました。

class RassoReg:

    def __init__(self, labmda_ = 1.0, max_iter = 1000):
        self.lambda_ = labmda_
        self.w_ = None
        self.max_iter = max_iter

    def thresholding(self, q, r,n):
        #wkが正の場合
        if q > self.lambda_ * n:
            return (q - self.lambda_*n) / r
        #wkが負の場合
        elif q < -self.lambda_ * n:
            return (q + self.lambda_*n) / r
        #その他
        else:
            return 0

    def _update(self, n, d, X, y):
        #バイアス(w0)の重みを更新
        self.w_[0] = (y - X[:,1:]@self.w_[1:]).sum() / n
        #バイアスの値がサンプルデータ数分だけ並んだベクトルを計算用に作成
        w0vec = np.ones(n) * self.w_[0]
        for k in range(d-1):
            ww = self.w_[1:]
            #バイアス更新対象の重みを計算式から省く(要数式参照)
            ww[k] = 0
            X_tmp = X[:,1:]
            #更新式の分子の一部
            q = (y - w0vec - X_tmp@ww)@X_tmp[:,k]
            #更新式の分母の一部
            r = X_tmp[:,k]@X_tmp[:,k]
            #重みの値を更新
            self.w_[k + 1] = self.thresholding(q, r,n)

    def fit(self, X, y):
        #バイアス計算用に1列目に1を挿入
        X = np.insert(X, 0, 1, axis=1)
        #nは行数（サンプルデータ数）、dは次元数（説明変数の数）
        n, d = X.shape
        #重みの初期値は0
        self.w_ = np.zeros(d)
        #指定回数繰り返し
        for _ in range(self.max_iter):
            self._update(n, d, X, y)

    def predict(self, X):
        return X@self.w

検証

ボストン住宅価格のデータセットを用いて検証を行います。ボストン住宅価格のデータセットの中身についてはこちらの記事(重回帰分析を理解する（自力実装編）で説明しているので、中身を知りたい方はご参照ください。

from sklearn.datasets import load_boston
import pandas as pd
from sklearn.preprocessing import StandardScaler

#データの読み込み
boston = load_boston()

#一旦、pandasのデータフレーム形式に変換
df = pd.DataFrame(boston.data, columns=boston.feature_names)

#目的変数（予想したい値）を取得
target = boston.target
df['target'] = target

X = df.iloc[:,:-1].values
X = StandardScaler().fit_transform(X)
y = df['target'].values

スクラッチ実装の結果

スクラッチ実装したラッソ回帰にボストン住宅価格のデータセットを当てはめるとこのようになりました。

linear = RassoReg()

linear.fit(X,y)
#切片
print(linear.w_[0])
#回帰係数
print(linear.w_[1:])

切片と回帰係数を出力します。

22.532806324110688
[ 0.          0.          0.          0.          0.          2.71310728
  0.          0.          0.          0.         -1.34349862  0.18079388
 -3.54361166]

sklearnの結果

sklearnのラッソ回帰にボストン住宅価格のデータセットを当てはめるとこのようになりました。

import sklearn.linear_model as lm
lasso = lm.Lasso(alpha=1.0, max_iter=1000, tol=0.0)
lasso.fit(X, y)
print(lasso.intercept_)
print(lasso.coef_)