PyTorchで線形回帰+L1/L2正則化

# はじめに
PyTorchを勉強しようと思ってUdemyで見つけたPyTorchの入門コースを受講しました。

**PyTorch Boot Camp : Python AI PyTorchで機械学習とデータ分析完全攻略**
https://www.udemy.com/course/python-pytorch-facebookai/

とてもよい感じだと思いました（小並感）。せっかく勉強したので、復習としてPyTorchを使って一番簡単な線形回帰をやってみたいと思います。
Pythonのバージョンは3.7.9、PyTorchは1.7.1です。

# 使用するデータセット・モデル
今回使用するのは以下のようなボストン住宅価格のデータセットです。
![スクリーンショット 2021-02-20 21.16.08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/6c65bfc6-09ca-9ee4-7e22-43862419eabc.png)

「CRIM」~「LSTAT」列を説明変数にして、目的変数である住宅価格「MEDV」の値を予測する線形回帰モデルを作成します。線形回帰なので以下のようなモデルで各説明変数に対する回帰係数$w_1, \ldots, w_{13}$と切片$b$の値を推定することで予測モデルを作成します。

```math
MEDV_i = b + w_1\times CRIM_i +  w_2\times ZN_i + \cdots + w_{13}\times LSTAT_i
```
添字$i$はデータの番号を表しています。

PyTorchでは「torch.nn.Linear」というのを使ってこのモデルを表現します。これはデータを線形変換する、ニューラルネットワークにおける全結合層を作るためのクラスです。上記の線形回帰モデルをニューラルネットワーク風に表現すると以下のような感じですね。

<img width="300" alt="nn.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/e8e82519-a301-8286-c248-fb73658117ad.png">

上の絵の通り、線形回帰モデルは入力層のノードが13、出力層のノードが1の2層のニューラルネットワークを作ることと同じになっています。中間層も活性化関数も無いのでニューラルネットワークと言って良いかは分かりませんが。今回は上の絵のようなモデルに対して、PyTorchを使って最適なパラメータ$w_1, \ldots, w_{13}$とバイアス$b$を学習することになります。

使用するデータセットはトレーニング用とテスト用で分けてcsvファイルで保存しておきます。

```python
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split

# データセット読み込み
boston = load_boston()
data = pd.DataFrame(boston.data, columns = boston.feature_names)
target = pd.DataFrame(boston.target, columns = ['MEDV'])
data = data.join(target)

# データ分割
train_data, test_data = train_test_split(data, test_size = 0.2, random_state = 1)

# データ保存
train_data.to_csv('boston_train.csv', index = False)
test_data.to_csv('boston_test.csv', index = False)
```

# PyTorchで実装

### モデル

```python
import torch.nn as nn

class LinearModel(nn.Module):
    def __init__(self, n_input, n_output):
        super().__init__()
        self.fc = nn.Linear(n_input, n_output, bias = True)
        
    def forward(self, x):
        predicted = self.fc(x)
        return predicted

model = LinearModel(n_input = 13, n_output = 1) # LinearModelのインスタンス作成
```
ここではnn.Moduleを継承し、予測に使うモデルのクラスLinearModelを定義しています。\_\_init_\_部分ではモデルの構造を作っています。n_inputで入力層のノードの数、n_outputで出力層のノードの数を指定し、nn.Linearでネットワークを作っています。複数の層を重ねる場合は、前層の出力ノードの数と次の層の入力ノードの数を揃える必要があります。biasはデフォルトでTrueなので省略可能です。forwardでは\_\_init_\_で定義したネットワークにデータxを入れて予測値の計算を行います。

最後にLinearModelのインスタンスを作成しています。入力層のノードの数を13、出力層のノードの数を1として引数で指定しています。この段階でモデルのパラメータの初期値が割り当てられているので、以下のようにLinearModelのインスタンスに特徴量を渡すことで予測値を算出することができます。

```python
features = torch.FloatTensor([1]*13) # 全ての値が1の特徴量を作成
model(features) # 予測値の計算
```
![スクリーンショット 2021-02-23 1.44.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/c8d724eb-3f92-cfae-4089-ffff0c0ad93f.png)

パラメータの値は以下のように確認できます。

```python
# 初期値確認
model.state_dict()
```
![スクリーンショット 2021-02-23 1.47.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/2b44ec28-2dc9-e469-35f0-c1502efb69f4.png)

### Dataset

```python
import pandas as pd
import torch
from torch.utils.data import Dataset
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

class BostonDataset(Dataset):
    def __init__(self):
        Boston = pd.read_csv('boston_train.csv')
        X = Boston.drop(columns = 'MEDV').values
        y = Boston['MEDV'].values
        
        scaler.fit(X)
        X = scaler.transform(X)
        
        self.X_train = torch.FloatTensor(X)
        self.y_train = torch.FloatTensor(y).view(-1, 1)
        
        self.data_length = len(self.y_train)
        
    def __getitem__(self, index):
        return self.X_train[index], self.y_train[index]
    
    def __len__(self):
        return self.data_length
```
Datasetを継承し、元データ（今回はboston_train.csv）から特徴量（説明変数）とそれに対応するラベル（目的変数）を1組づつ取り出すことのできるBostonDatasetクラスを作成します。\_\_init_\_でcsvファイルを開いて、特徴量とラベルをPyTorchで使用するTensor型に変換しています。また、説明変数の標準化も行なっています。Datasetでは以下のように\_\_getitem_\_と\_\_len_\_を実装する必要があります。

**\_\_len\_\_**：データの長さ（今回はcsvファイルの行数）を返すようにする。
**\_\_getitem\_\_**：インデックス番号に対応する特徴量とそのラベルを返すようにする。

このように実装することで、object[index]で対応する番号のデータにアクセスできるようになり、len(object)でデータの長さをを取得できるようになります。Datasetは後述のDataLoaderを作るときに使用します。

### DataLoader

```python
from torch.utils.data import DataLoader

batch_size = 64  # バッチサイズ
data_set = BostonDataset() #BostonDatasetのインスタンス作成
trainloader = DataLoader(dataset = data_set, batch_size = batch_size, shuffle = True) #DataLoaderを作成
```
DataLoaderを使用することでDatasetからバッチサイズ分のデータを取り出すことができるようになります。DataLoader()の引数に先ほど作成したBostonDatasetのインスタンス、バッチサイズを指定することで作成します。また、shuffle = Trueを指定することでデータをシャッフルしてくれます。

作成したDataLoaderの挙動を確認してみます。

```python
print('Data length:',len(data_set))
for i, (X, y) in enumerate(trainloader):
    print(f'Batch: {i}, X size: {X.size()}, label size: {y.size()}')
```

![スクリーンショット 2021-02-23 1.21.55.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/4cb40dce-d0f2-0d8a-fc87-461a1581c0af.png)
上記のようにDataLoaderをfor分で処理すると、バッチサイズ分の特徴量とそのラベルを順番に取り出すことができます。総データ数が404、バッチサイズを64にしているので、イテレーション数は7で最後はデータ数が20になっています。

### 損失関数、最適化アルゴリズム

```python
loss_function = nn.MSELoss()
optimizer = torch.optim.Adam(params = model.parameters(), lr = 0.01)
```
損失関数と最適化アルゴリズムを指定します。今回は損失関数として平均二乗誤差、最適化アルゴリズムはAdamを用いることにします。lrは学習率でここでは0.01としておきます。

### 学習の実行
準備が整ったので学習を実行します。今回はエポックを800に設定してみます。設定したエポック数とDataLoaderに対してfor文を実行して学習・パラメータの更新を行います。

```python
epochs = 800 #エポックを800に設定
loss_list = [] #損失関数の値を保存するためのリスト

for epoch in range(epochs):
    for X, y in trainloader:
        
        y_pred = model(X) # 予測値の計算
        loss = loss_function(y_pred, y) # 損失関数の値を計算
        
        optimizer.zero_grad() # 勾配を初期化
        loss.backward() # 勾配を計算
        optimizer.step() # パラメータを更新
        
        loss_list.append(loss.detach()) # 損失関数の値を保存
```
PyTorchでの学習は上記のコードのように「予測値の計算→損失関数の値を計算→勾配初期化→勾配計算→パラメータ更新」という順番で実行していきます。モデルが変わっても基本的には同様の流れで学習を実行することができます。

### 学習結果の確認
損失関数の値の推移、最終的な損失関数の値を確認しておくことにします。

```python
import matplotlib.pyplot as plt

print(loss_list[-1])
plt.plot(loss_list)
plt.xlabel('Training Iteration') 
plt.ylabel('Loss') 
plt.show()
```

![スクリーンショット 2021-02-26 22.11.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/464d9dcc-b56b-6aed-64b2-059adb1834a8.png)



損失関数の値の推移を見ると学習はほぼ収束していることがわかります。学習終了時点での損失関数の値(MSE)は約13.7になっていました。

### モデルの評価
あらかじめ分けておいたテスト用のデータを使って、テスト用データに対するMSEを計算してモデルの性能を評価することにします。

```python
# データ読み込み
test_data = pd.read_csv('boston_test.csv')

X_test = test_data.drop(columns = 'MEDV').values
X_test = scaler.transform(X_test)
X_test = torch.FloatTensor(X_test)

y_test = test_data['MEDV'].values
y_test = torch.FloatTensor(y_test).view(-1, 1)

# テスト用データでMSEを計算
model.eval()
with torch.no_grad():
    y_pred = model(X_test)
    loss = loss_function(y_pred, y_test)
    print(loss)
```
![スクリーンショット 2021-02-25 22.53.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/687f9782-546a-6471-2242-4fa7d7a9ef30.png)

モデルの評価を行うときはmodel.eval()で評価モードに切り替えています。また、評価時には勾配の計算は必要ないので、とりあえずtorch.no_grad()をするようにしています。今回のモデルのテスト用データに対するMSEは約23.2でした。

# L1/L2正則化
ついでにL1正則化とL2正則化をPyTorchで試してみることにします。正則化は損失関数に以下の①L1ノルムもしくは②L2ノルムの二乗を加えれば良いので、学習を行うときのコードを書き換えるだけでOKです。$\alpha$は正則化パラメータです。

```math
①:\alpha\sum_{j}|\beta_j|
```

```math
②:\alpha\sum_{j}\beta_j^2
```
パラメータ学習時のコード以外はこれまでと同じものを利用しています。正則化するときは目的変数も標準化した方が良いかもしれないですが、今回はここまでと同じで標準化は説明変数のみでいきます。また、以下のコードでパラメータを初期化してから学習を行いました。

```python
# パラメータ初期化
model.fc.reset_parameters()
```

### L1正則化（Lasso回帰）
L1正則化はパラメータのL1ノルムを損失関数に加えれば良いので、これまで同様にMSEを計算してそこにL1ノルムを足す、というコードを書き加えています。また、torch.norm()を使ってノルムの計算行なっています。スパースなモデルが作れているか確認するために、正則化パラメータ$\alpha$は1に設定して強めに正則化を行います。

```python
epochs = 800
loss_list = []
alpha = 1 # 正則化パラメータ

for epoch in range(epochs):
    for X, y in trainloader:
        
        y_pred = model(X)
        loss = loss_function(y_pred, y)
        
        # パラメータのL1ノルムを損失関数に足す
        l1 = torch.tensor(0., requires_grad=True)
        for w in model.parameters():
            l1 = l1 + torch.norm(w, 1)
        loss = loss + alpha*l1
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        loss_list.append(loss.detach())
```
損失関数の推移とトレーニングに用いたデータに対する損失関数のの値を確認します。ここでは損失関数=MSE+L1ノルムの値です。↓

![スクリーンショット 2021-02-27 23.41.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/9d814fed-0437-8826-f342-5ff4f0297fb7.png)

学習は収束していることがわかります。推定したパラメータの値をプロットして確認してみます。

```python
# パラメータの値を取得(バイアスは無し)
weight = model.state_dict()['fc.weight']

plt.scatter(list(range(1,14)), weight[0])
plt.axhline(0, color = "magenta")
plt.xticks(list(range(1,14)))
plt.xlabel('Parameter index', fontsize=14) 
plt.ylabel('Parameter magnitude', fontsize=14) 
plt.show()
```
![スクリーンショット 2021-02-27 23.56.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/ff946e7a-f317-15ae-63d8-89c36310d785.png)

縦軸に学習したパラメータの値をとって横にパラメータを並べています。いくつかのパラメータは0になっていることがわかります。L1正則化はこんな感じで行えます。

### L2正則化（Ridge回帰）
L2正則化は正則化項としてL2ノルムの二乗を加えれば良いので、L1正則化と同様に学習時のコードを以下のように書き換えればOKです。正則化パラメータ$\alpha$は0.01にしておきます。

```python
epochs = 800
loss_list = []
alpha = 0.01 # 正則化パラメータ

for epoch in range(epochs):
    for X, y in trainloader:
        
        y_pred = model(X)
        loss = loss_function(y_pred, y)
        
        # パラメータのL2ノルムの二乗を損失関数に足す
        l2 = torch.tensor(0., requires_grad=True)
        for w in model.parameters():
            l2 = l2 + torch.norm(w)**2
        loss = loss + alpha*l2
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        loss_list.append(loss.detach())
```

損失関数の推移とトレーニングに用いたデータに対する損失関数の値を確認します。ここでは損失関数=MSE+L2ノルムの二乗の値です。↓

![スクリーンショット 2021-02-26 22.40.23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/25e0b84f-1753-42f3-af7f-59357e5ecd91.png)

また、テスト用データに対するMSEの値も同様に確認します。↓

![スクリーンショット 2021-02-26 22.42.22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697043/df278c6f-e25e-30e1-cd69-c2501369b328.png)

正則化を行なっていないモデルと比べてテスト用データに対するMSEの値はちょっとだけ小さくなりました。正則化の効果があったとは言えないと思いますが、とりあえずL2正則化はこんな感じで行うことができます。また、L2正則化は最適化アルゴリズムを設定するときにtorch.optim.Adam(weight_decay=0.01)のように、weight_decayに正則化パラメータの値を入れることでも行えるようです。

# 終わりに
PyTorchの使い方が掴めてきたので、さらに複雑なモデルとかCNNとかもやっていきたいです。あと、どうやって勾配の計算してるかとかちゃんと理解できてない部分も結構あるので勉強していきたいです。

# 参考文献
・公式ドキュメント
https://pytorch.org/docs/stable/index.html

・正則化のやり方とか
https://stackoverflow.com/questions/42704283/adding-l1-l2-regularization-in-pytorch