#まえがき
ラッソ回帰についてまとめます。一応前回までの記事を読んでいる前提で書いています。

ラッソ回帰

ラッソ回帰において最小化する二乗誤差は以下の式で表されます。

E(\mathbf{w}) = ||\mathbf{y} - \mathbf{\tilde{X}\mathbf{w}}||^2 + \lambda ||w ||_1

ここで$|w|_1$は正則化項で、$|w|_1= \sum ^d _{i=1} |w _{i}|$です。
要するにwの各値の合計。二乗ではなく1乗で合計する。プラマイゼロになっちゃ困るので、それぞれ絶対値をとっています。
#ラッソ回帰を計算する上での致命的な問題
ラッソ回帰を解くためのアルゴリズムはいくつかありますが、今回は座標降下法で計算してみましょう。
簡単に言えば

変数を一つ選んで偏微分係数を0になる値に合わせる。
さらにもう一つ変数を選び、偏微分係数が0になる値に合わせる。
同じ工程を全変数に対して行う。

そこまでを一周として、学習が収束するまで何周も繰り返す方法です。

さて、この式の問題点は、今回の$E(w)$は「微分不可能である」という点です。絶対値が入っていると関数が連続にならないためです。ならば場合分けしなければいけない。

数式を書いてみましょう

まず、$w_0$で微分してみましょう。これは絶対値が含まれていないため、普通に偏微分して$w_0$を求めることができます。

\frac{\partial E}{\partial w_0} = -\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})

これが=0になるので、$w_0$について整理して

w_0 = {1\over n}\sum^n _ {i = 1}(y_i - \sum^d _ {j = 1}x_{ij}w_{j})

と求まりました。
$w_1$以降が厄介です。絶対値の微分を確認しましょう。

d^+(|x|) = 1\\
d^-(|x|) =-1

このように0の左側と右側で分かれます。
同様にE(w)を$w_k (k \neq 0)$で微分すると以下のように分かれます。

\frac{\partial E}{\partial w_k} =
\begin{cases}
 -\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} + \lambda & (w_k > 0)\\
 -\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} - \lambda
& (w_k < 0)
\end{cases}

どちらも=0として整理すると

w_k = \begin{cases}
{{\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} - \lambda} \over {\sum^n _{i = 1}x^2_{ik}w_k}} & (w_k > 0)\\
{{\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} + \lambda} \over {\sum^n _{i = 1}x^2_{ik}w_k}} & (w_k < 0)

\end{cases}

と求まりました。
ただし、これは$w_k$の条件を正しく満たす必要があります。上の式の条件を満たさない場合は$w_k = 0$となります。
書き直してみます。

w_k = \begin{cases}
{{\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} - \lambda} \over {\sum^n _{i = 1}x^2_{ik}w_k}} & ({\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik}  } > \lambda)\\
0 & (\lambda < {\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik}  } < \lambda)\\
{{\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} + \lambda} \over {\sum^n _{i = 1}x^2_{ik}w_k}} & (\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} < - \lambda)

\end{cases}

$w_k = 0$というのは、その要素は計算に含めないほうがいい、という意味です。実際に計算を行うとほとんどの要素が0になったりするので疎行列の扱いを覚えた方がいいようです。
後で計算を楽にするために以下のような工夫をしましょう。

w_k = {{S({\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik}} , \lambda) } \over {\sum^n _{i = 1}x^2_{ik}w_k}}

という関数$S$に置き換えて（ソフト閾値関数といいます）表すと、

{S({\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik}} , \lambda) }\\ = \begin{cases}
{{\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} - \lambda} } & ({\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik}  } > \lambda)\\
0 & (-\lambda < ({\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik}  } < \lambda )\\
{{\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} + \lambda} } & (\sum^n_{i=1}(y_i - w_0 - \sum^d_{j = 1}x_{ij}w_{j})x_{ik} < - \lambda)
\end{cases}

という風に整理できました。このSをsoft_thresholdingという名前で実装します。

実装してみよう

今回サンプルとして扱うデータは以下のサイトから利用します。
https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/
ここから winequality-red.csvというファイルをダウンロードして、同じディレクトリに入れました。このデータはワインの各種ステータスとその品質を比較することができます。

まずこのcsvファイルを読み込んでみましょう。
以下のように記述します。

import numpy as np
import csv

Xy = []
with open("winequality-red.csv") as fp:
    for row in csv.reader(fp, delimiter=";"):
        Xy.append(row)
Xy = np.array(Xy[1:], dtype=np.float64)

こういうデータが

こういうデータになりました。

訓練用のデータとテスト用のデータに分けます。

np.random.seed(0)
np.random.shuffle(Xy)
train_X = Xy[:-1000, :-1] # 1件を除いて訓練用のデータに格納する。
train_y = Xy[:-1000, -1] # 訓練データの品質一覧。
test_X = Xy[-1000, :-1] # テストデータ(1件)のパラメータ。これを使って品質を計算する。
test_y = Xy[-1000, -1] # テストデータの品質(1件)。これを当てられたらOK。

あとは処理を入力します。

def soft_thresholding (x, y):
    return np.sign(x) * max(abs(x) - y, 0)

class Lasso:
    def __init__(self, lambda_, tol = 0.0001, max_iter = 10000):
        self.lambda_ = lambda_
        self.tol = tol
        self.max_iter = max_iter
        self.w_ = None


    def fit(self, X, t):
        n, d = X.shape
        self.w_ = np.zeros(d + 1)
        avgl1 = 0.
        for _ in range(self.max_iter): # max_iter = 10000回繰り返す。
            avgl1_prev = avgl1
            self._update(n, d, X, t)
            avgl1 = np.abs(self.w_).sum() / self.w_.shape[0]
            if abs(avgl1 - avgl1_prev) <= self.tol:
                break # 誤差がtol = 0.0001以下になったら終了。

    def _update(self, n, d, X, t):
        self.w_[0] = (t - np.dot(X, self.w_[1:])).sum() / n
        w0vec = np.ones(n) * self.w_[0]
        for k in range(d):
            ww = self.w_[1:]
            ww[k] = 0
            q = np.dot(t - w0vec - np.dot(X, ww), X[:, k])
            r = np.dot(X[:, k], X[:, k])
            self.w_[k + 1] = soft_thresholding(q / r, self.lambda_)

    def predict(self, X):
        if X.ndim == 1:
            X = X.reshape(1, X.shape[0])
            Xtil = np.c_[np.ones(X.shape[0]), X]
        return np.dot(Xtil, self.w_)

for lambda_ in [1., 0.1, 0.01]: # 三種類のハイパーパラメータを試す。
    model = Lasso(lambda_)
    model.fit(train_X, train_y)
    y = model.predict(test_X)
    print("--- ハイパーパラメータlambda = {} ---".format(lambda_))
    print("係数w:")
    print(model.w_)
    mse = ((y - test_y)**2).mean()
    print("予想値: {}".format(y[0]))
    print("実際の値:{}".format(test_y))
    print("誤差E(w): {:.3f}".format(mse))

出力は以下のようになります。

--- ハイパーパラメータlambda = 1.0 ---
係数w:
[ 5.6360601  0.        -0.         0.         0.        -0.
 -0.        -0.        -0.        -0.         0.         0.       ]
予想値: 5.636060100166945
実際の値:7.0
誤差E(w): 1.860
--- ハイパーパラメータlambda = 0.1 ---
係数w:
[ 6.13255766  0.         -0.75763602  0.27031068  0.         -1.97134698
 -0.         -0.          0.          0.          0.          0.        ]
予想値: 5.844448204297558
実際の値:7.0
誤差E(w): 1.335
--- ハイパーパラメータlambda = 0.01 ---
係数w:
[ 5.76652084 -0.         -1.51392331 -0.03358247  0.         -3.10993286
 -0.         -0.          0.21500579  0.          1.08349525  0.        ]
予想値: 6.165807417800908
実際の値:7.0
誤差E(w): 0.696

このように計算できました。
ハイパーパラメータが高いほど$w_k$が「ないほうがマシ」と判定されて疎行列になっていくこと、逆にハイパーパラメータが低いほど過学習のリスクが高まることがわかります。

まとめ

ラッソ回帰は絶対値が含まれて計算が厄介な数式でしたが、数学の力で特別なライブラリなど使用せずに実装できました。
ハイパーパラメータの設定によって結果が変わるので注意が必要です。

『機械学習のエッセンス』という本を参考にしました。この本ではラッソ回帰のコードには誤植があるので気を付けてください。
前回:多項式回帰、過学習
https://qiita.com/NNNiNiNNN/items/d87990a6eef72a3815a3