線形回帰①

# 線形回帰とは
### 線形回帰モデル
線形回帰モデルの一般的な形式は以下のように表されます：
$$\hat{y} = \theta_0 + \theta_1 x_1 + \cdots + \theta_n x_n$$
- $\hat{y}$:応答変数（または従属変数、目的変数）と呼ばれる、予測したい目標の値。
- $x_i$:予測変数（または独立変数、特徴量）呼ばれる、応答変数を予測するのに使う情報。
- $n$:予測変数の数
- $\theta_j$:モデルにおける$j$番目のパラメータ。$j$番目の予測変数に対応する。

これをベクトル形式で書き表すと以下のようになります：
$$ \hat{y} = \boldsymbol{\theta}^\mathsf{T}\boldsymbol{x} $$

$\boldsymbol{\theta}$と$\boldsymbol{x}$はパラメータと予測変数を列ベクトルにしたもので、これらのベクトルから転置と内積を使って線形回帰モデルを表すことができます。
※バイアス項に対応する$x_0$は常に1

線形回帰モデルの目的は、訓練データを使用してこれらのパラメタ$\boldsymbol{\theta}$を学習し、予測変数$\hat{y}$と実際の値$y$との間の誤差を最小化することです。この誤差の測定方法は様々なものがありますが、代表的なものが以下で表される平均二乗誤差（Mean Square Error:MSE）です。
$$MSE=\frac{1}{m}\sum_{i=1}^{m}\left( \boldsymbol{\theta}^\mathsf{T}\boldsymbol{x}^{(i)} - y^{(i)} \right)^2$$
$\boldsymbol{x}$と$y$の右肩の添字$(i)$は$i$番目のデータを表し、全体では$m$個のデータがあることを意味しています。すなわち最初の式では個々のデータ点での予測モデルを意味していましたが、MSEを考える際は与えられた訓練データについての予測誤差の二乗を合計しデータ個数で割ったもの（平均）を計算します。

この$MSE$を最小化するようなパラメータ$\boldsymbol{\theta}$を求めたいのですが、求め方として代表的なものが最小二乗法と勾配降下法です。

### 最小二乗法
$MSE$を最小化するために、$MSE$をパラメタ$\boldsymbol{\theta}$で微分して0となる値を求めるのが最小二乗法ですが、まずは$MSE$を微分しやすいように変形していきます。

線形回帰をベクトルで表した式を$m$個のデータセットについて列ベクトルで並べると
```math
\begin{align}
\begin{pmatrix}
 \hat{y}^{(1)} \\ \hat{y}^{(2)} \\ \vdots \\ \hat{y}^{(m)} 
\end{pmatrix}
&=
\begin{pmatrix}
 \boldsymbol{\theta}^\mathsf{T}\boldsymbol{x}^{(1)} \\ 
 \boldsymbol{\theta}^\mathsf{T}\boldsymbol{x}^{(2)} \\ 
 \vdots \\
 \boldsymbol{\theta}^\mathsf{T}\boldsymbol{x}^{(m)} 
\end{pmatrix} \\
&=
\begin{pmatrix}
 \boldsymbol{x}^{(1)\mathsf{T}} \\ \boldsymbol{x}^{(2)\mathsf{T}} \\ \vdots \\ \boldsymbol{x}^{(m)\mathsf{T}} 
\end{pmatrix}\boldsymbol{\theta} \\
&=
\begin{pmatrix}
 x_{0}^{(1)} & x_{1}^{(1)} & \cdots & x_{n}^{(1)} \\ 
 x_{0}^{(2)} & x_{1}^{(2)} & \cdots & x_{n}^{(2)} \\ 
 & \vdots \\ 
 x_{0}^{(m)} & x_{1}^{(m)} & \cdots & x_{n}^{(m)}
\end{pmatrix}\boldsymbol{\theta}
\end{align}
```
左辺の$m$個の予測値のベクトルを$\hat{\boldsymbol{y}}$, 右辺の予測変数の行列を$\boldsymbol{X}$で表すと、
```math
\begin{align}
\hat{\boldsymbol{y}} = \boldsymbol{X}\boldsymbol{\theta}
\end{align}
```
となります。これを使って$MSE$の式を変形していくと、
```math
\begin{align}
MSE
&= \frac{1}{m}\sum_{i=1}^{m}\left( \boldsymbol{\theta}^\mathsf{T}\boldsymbol{x}^{(i)} - y^{(i)} \right)^2 \\
&= \frac{1}{m}\left( \boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y} \right)^\mathsf{T}\left( \boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y} \right) \\
&= \frac{1}{m}\left( \left(\boldsymbol{X}\boldsymbol{\theta}\right)^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - \left(\boldsymbol{X}\boldsymbol{\theta}\right)^\mathsf{T}\boldsymbol{y}
 - \boldsymbol{y}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 + \boldsymbol{y}^\mathsf{T}\boldsymbol{y} \right)\\
&= \frac{1}{m}\left( \left(\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\right)\boldsymbol{X}\boldsymbol{\theta}
 - \left(\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\right)\boldsymbol{y}
 - \boldsymbol{y}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 + \boldsymbol{y}^\mathsf{T}\boldsymbol{y} \right)\\
&= \frac{1}{m}\left( \boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - \boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 - \boldsymbol{y}^\mathsf{T}\left(\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\right)^\mathsf{T}
 + \boldsymbol{y}^\mathsf{T}\boldsymbol{y} \right)\\
&= \frac{1}{m}\left( \boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - \boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 - \left(\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}\boldsymbol{y}\right)^\mathsf{T}
 + \boldsymbol{y}^\mathsf{T}\boldsymbol{y} \right)\\
&= \frac{1}{m}\left( \boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - 2\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 + \boldsymbol{y}^\mathsf{T}\boldsymbol{y} \right)
\end{align}
```
と変形できるため、これを$\boldsymbol{\theta}$で微分して
```math
\begin{align}
\frac{\partial }{\partial \boldsymbol{\theta}} MSE 
&= \frac{1}{m} \frac{\partial }{\partial \boldsymbol{\theta}}\left(
\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - 2\boldsymbol{\theta}^\mathsf{T}\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 + \boldsymbol{y}^\mathsf{T}\boldsymbol{y}
 \right) \\
&= \frac{1}{m} \left(\left(
\boldsymbol{X}^\mathsf{T}\boldsymbol{X}+\left(\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\right)^\mathsf{T}\right)\boldsymbol{\theta}
 - 2\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 \right) \\
&= \frac{1}{m} \left(
2\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - 2\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 \right)
\end{align}
```
※ベクトルの微分の公式$\frac{\partial }{\partial \boldsymbol{A}}\boldsymbol{A}^\mathsf{T}\boldsymbol{B}=\boldsymbol{B}$および$\frac{\partial }{\partial \boldsymbol{A}}\boldsymbol{A}^\mathsf{T}\boldsymbol{B}\boldsymbol{A}=\left(\boldsymbol{B}+\boldsymbol{B}^\mathsf{T}\right)\boldsymbol{A}$を使っています。

$MSE$を微分して0となる$\boldsymbol{\theta}$を求めると
```math
\begin{align}
\frac{1}{m} \left(
2\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - 2\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 \right)
= 0 \\
\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - \boldsymbol{X}^\mathsf{T}\boldsymbol{y}
= 0 \\
\boldsymbol{\theta}
= \left( \boldsymbol{X}^\mathsf{T}\boldsymbol{X} \right)^{-1}\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
\end{align}
```
これが最小二乗法により求められた、$MSE$を最小にするパラメタ$\boldsymbol{\theta}$です。




### 勾配降下法
勾配降下法でも、$MSE$を最小にする$\boldsymbol{\theta}$を微分により求めますが、そのアプローチは異なります。まず適当なパラメタの初期値$\boldsymbol{\theta_0}$を決めます。その後$MSE$の勾配$\nabla_{\theta}MSE$を使ってパラメタを更新します。このステップを繰り返すことで最終的にパラメタ$\boldsymbol{\theta}$を求めることができます。
```math
\begin{align}
\boldsymbol{\theta^{next}} = \boldsymbol{\theta}-\eta\nabla_{\theta}MSE \\
\nabla_{\theta}MSE = \frac{1}{m} \left(
2\boldsymbol{X}^\mathsf{T}\boldsymbol{X}\boldsymbol{\theta}
 - 2\boldsymbol{X}^\mathsf{T}\boldsymbol{y}
 \right)
\end{align}
```
$\eta$は学習率と呼ばれる値で各ステップでの$\theta$を更新する大きさを決めます。$\eta$が大きければ$\theta$は大きく更新されていくため学習は速くなりますが、学習が収束しない恐れがあります。$\eta$を小さくすると、最終的には適切な解に収束する可能性が高いですが、学習が終わるまでに大きな時間がかかってしまいます。学習率はパラメタ$\theta$とは別に設定する必要があるため、いくつかの学習率で学習させて結果を比較することが一般的です。

$\nabla_{\theta}MSE$に含まれる$\boldsymbol{X}$はn個の成分をもつ予測変数をデータセットm個分だけ並べたものですが、このように全データセットを使って勾配を計算しパラメタを更新していく方法をバッチ勾配降下法と呼びます。バッチ勾配降下法ではデータ数が大きくなるにつれ計算コストが増大してしまうという欠点があります。一方同じ勾配法でも、確率的勾配降下法と呼ばれる方法では、データセット$\boldsymbol{X}$から一つだけデータ$\boldsymbol{x}^{(i)}$を無作為に選び勾配を計算します。そして一つのデータから計算した勾配でパラメータを更新し、また別のデータセットをランダムに選択し勾配計算とパラメタ更新を繰り返します。m個のデータセットに対し、重複がないように各データが1回ずつ選択されるようにm回の計算を行い、これを1エポック（epoch）と呼びます。計算はパラメータ更新が収束するまで数十エポック程度行います。確率的勾配降下法は計算効率が良いという長所がありますが、1データの勾配のみからパラメタを更新するので収束が不安定になる可能性がある欠点があります。この対処法として、学習率を徐々に小さくしていく方法があります。

バッチ勾配降下法と確率的勾配降下法の長所を組み合わせた方法がミニバッチ勾配降下法です。ミニバッチ勾配降下法では全体のデータセットからミニバッチと呼ばれるいくつかのデータを選択し勾配の計算とパラメータの更新を行います。バッチ勾配降下法と比較し計算効率がよく、確率的勾配降下法よりも収束が安定する特徴があります。

# 正則化
最小二乗法でも勾配法でも平均二乗誤差を最小化することを目指しますが、学習データに対する誤差が小さくなりすぎると学習に使っていないデータに対する予測誤差が大きくなってしまう事があり、これを過学習と呼びます。実際の機械学習が役に立つ場面では与えられたデータを説明するだけでなく未知のデータ、新しいデータに対する予測も重要になってくるため過学習は解決すべき課題と言えます。正則化とは過学習を防ぐための手法であり、モデルのパラメタに制約を設けます。

### L1正則化
L1正則化（線形回帰においてはLasso回帰）と呼ばれる手法は、最小二乗誤差にパラメタの絶対値の和を加えたものをコスト関数と定義しこの勾配を用いてパラメタを更新する手法です。具体的なコスト関数$J(\theta)$は
```math
\begin{align}
J(\boldsymbol{\theta}) = MSE + \alpha\sum_{i=1}^n|\theta_i|
\end{align}
```
と定義されます。追加した$\alpha\sum_{i=1}^n|\theta_i|$は正則化項と呼ばれ、$\alpha$はハイパーパラメタで正則化の強さを決定します。L1正則化のコスト関数は絶対値を含むためそのままでは勾配が計算できませんが、$\rm{sign(\theta)}$を用いて下記のように計算することで勾配降下法は機能するようになります。
```math
\begin{align}
\frac{\partial }{\partial \boldsymbol{\theta}}J(\boldsymbol{\theta}) = \nabla_{\theta}MSE + \alpha
\begin{pmatrix}
 \rm{sign(\theta_1)} \\ \rm{sign(\theta_2)} \\ \vdots \\ \rm{sign(\theta_n)} 
\end{pmatrix} \\
\rm{sign(\theta_i)} = 
   \begin{cases}
      -1  \;\; (\theta_i<0) \\
      0  \;\;(\theta_i=0) \\
      +1 \;\;(\theta_i>0)
   \end{cases}
\end{align}
```

### L2正則化
L2正則化（線形回帰においてはRidge回帰）もL1正則化と同様にコスト関数に正則化項を追加しますが、正則化項はパラメタの二乗和となります。
```math
\begin{align}
J(\boldsymbol{\theta}) = MSE + \frac{\alpha}{2}\sum_{i=1}^n \theta_i^2
\end{align}
```
L2正則化ではコスト関数がそのまま微分できるため勾配は下記のようになります。
```math
\begin{align}
\frac{\partial }{\partial \boldsymbol{\theta}}J(\boldsymbol{\theta}) = \nabla_{\theta}MSE + \alpha\boldsymbol{\theta}
\end{align}
```
L1正則化もL2正則化もどちらもデータに過学習しないようにペナルティを与えますが、L2正則化がどのパラメタにも平等にペナルティを課すのに対し、L1正則化は重要度が高い特徴量を残し重要度が低い特徴量を無視するようにペナルティを課します。同じ正則化でもこのような特徴の違いを理解することでデータセットの特徴からモデルを使い分けることができます。