#まえがき
ラッソ回帰についてまとめます。一応前回までの記事を読んでいる前提で書いています。
ラッソ回帰
ラッソ回帰において最小化する二乗誤差は以下の式で表されます。
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