線形回帰・リッジ回帰・ラッソ回帰 -計算方法-

線形回帰・リッジ回帰・ラッソ回帰を行う際の計算方法

* 参考にしたもの：[機械学習のエッセンス](https://www.sbcr.jp/product/4797393965/)の第05章「機械学習アルゴリズム」
* ざっくりとした考え方は[コッチ](https://qiita.com/kazuya_minakuchi/items/186b20efd0eefdeb026b)

# 線形回帰

yを予測するような$f(\mathbf{x})$を↓のように表し、できるだけ正確に予測できるよう$\mathbf{w}^t=(w_0, w_1, w_2, \cdots, w_d)$を求める

```math
f(\mathbf{x}) = w_0 + w_1x_1 + w_2x_2 + \cdots + w_d x_d = w_0 + \sum_{i=1}^{d} w_i x_i = \mathbf{w}^t \mathbf{x} 
```

* $y$：目的変数。真の値。コレを予測したい。
* $\mathbf{x}^t = (1, x_1, x_2, \cdots, x_d)$：説明変数。コレを使って$y$を予測したい。(最初の1は説明変数ではない。行列計算するために追加したもの)
  * この場合1つの説明変数はd次元
* $f(\mathbf{x})$：予測値

↓こういうの
![線形回帰.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/563901/09cd4916-4c69-1c03-abce-a7f3433b9a8d.png)

直線を引いたときに「予測と真の値の平均二乗誤差」

```math
E = \sum (y - f(x))^2
```

が最小になるよう、パラメータ(係数)を求める。
ここで、教師データ($\mathbf{x}$、$y$の組み合わせ)がn個あり、

```math
\mathbf{y} = 
\begin{pmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n
\end{pmatrix}
```

```math
\mathbf{X} = 
\begin{pmatrix}
{x_1}^t \\
{x_2}^t \\
\vdots \\
{x_n}^t
\end{pmatrix}
```

とすると。$E$は

```math
\begin{align}
E(\mathbf{w}) &= \sum_{i=1}^{n}{(y_i - f(\mathbf{x_i}))^2} \\
&= \sum_{i=1}^{n}{(y_i - \mathbf{w}^t \mathbf{x_i})^2} \\
&= ||\mathbf{y} - \mathbf{X} \mathbf{w}||^2 \\
&= (\mathbf{y} - \mathbf{X} \mathbf{w})^t (\mathbf{y} - \mathbf{X} \mathbf{w}) \\
&= \mathbf{y}^t \mathbf{y} - \mathbf{w}^t \mathbf{X}^t \mathbf{y} - \mathbf{y}^t \mathbf{X} \mathbf{w} + \mathbf{w}^t \mathbf{X}^t \mathbf{X} \mathbf{w} \\
\end{align}
```

と表せる。コレを最小にするような$\mathbf{w}$は、Eの勾配

```math
\nabla E(\mathbf{w}) = -2 \mathbf{X}^t \mathbf{y} + 2 \mathbf{X}^t \mathbf{X} \mathbf{w}
```

が0になる$\mathbf{w}$なので

```math
\begin{align}
-2 \mathbf{X}^t \mathbf{y} + 2 \mathbf{X}^t \mathbf{X} \mathbf{w} &= 0 \\
\mathbf{w} &= (\mathbf{X}^t \mathbf{X})^{-1} \mathbf{X}^t \mathbf{y}
\end{align}
```

となる。

# リッジ回帰

「(予測と真の値の平均二乗誤差)+<font color="Red">(ハイパーパラメータ)×(パラメータのL2ノルム)</font>」

```math
E(\mathbf{w}) = \sum_{i=1}^{n}{(y_i - f(\mathbf{x_i}))^2} + \lambda \sum_{i=0}^{d} w_i^2
```

が最小になるよう、パラメータ(係数)を求める(=L2正則化)
このとき

```math
\begin{align}
E(\mathbf{w}) &= ||\mathbf{y} - \mathbf{X} \mathbf{w}||^2 + \lambda ||\mathbf{w}||^2 \\
&= \mathbf{y}^t \mathbf{y} - \mathbf{w}^t \mathbf{X}^t \mathbf{y} - \mathbf{y}^t \mathbf{X} \mathbf{w} + \mathbf{w}^t \mathbf{X}^t \mathbf{X} \mathbf{w} + \lambda ||\mathbf{w}||^2 \\
\end{align}
```

となり、コレを最小にするような$\mathbf{w}$は、Eの勾配

```math
\nabla E(\mathbf{w}) = -2 \mathbf{X}^t \mathbf{y} + 2 \mathbf{X}^t \mathbf{X} \mathbf{w} + 2 \lambda \mathbf{w}
```

が0になる$\mathbf{w}$なので

```math
\begin{align}
-2 \mathbf{X}^t \mathbf{y} + 2 \mathbf{X}^t \mathbf{X} \mathbf{w} + 2 \lambda \mathbf{w} &= 0 \\
(\mathbf{X}^t \mathbf{X} + \lambda \mathbf{I}) \mathbf{w} &= \mathbf{X}^t \mathbf{y} \\
\mathbf{w} &= (\mathbf{X}^t \mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X}^t \mathbf{y}
\end{align}
```

となる。

# ラッソ回帰

「(予測と真の値の平均二乗誤差)+<font color="Red">(ハイパーパラメータ×パラメータのL1ノルム)</font>」が最小になるよう、パラメータ(係数)を求める。=L1正則化

```math
\begin{align}
E(\mathbf{w}) &= \frac{1}{2} \sum_{i=1}^{n}{(y_i - f(\mathbf{x_i}))^2} + \lambda \sum_{i=1}^{d} |w_i| \\
&= \frac{1}{2} \sum_{i=1}^{n}{(y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij})^2} + \lambda \sum_{i=1}^{d} |w_i| \\
&= \frac{1}{2} ||\mathbf{y} - \mathbf{X} \mathbf{w}||^2 + \lambda |\mathbf{w}|_1 \\
&= \mathbf{y}^t \mathbf{y} - \mathbf{w}^t \mathbf{X}^t \mathbf{y} - \mathbf{y}^t \mathbf{X} \mathbf{w} + \mathbf{w}^t \mathbf{X}^t \mathbf{X} \mathbf{w} + \lambda |\mathbf{w}|_1 \\
\end{align}
```

とおく。ラッソ回帰は$E(\mathbf{w})=0$となる解を手計算で求められないので、数値計算を行う。

## 数値計算(座標降下法:Coordinate Descent)

- やりたいこと
    - $E(w_0, w_1, w_2, \cdots, w_d)$を最小にする$w_0, w_1, w_2, \cdots, w_d$を求める
- 手順
    1. まず$w_0$を求める
    2. 次に、適当に$\mathbf{w}^{(0)} = (w_0, w_1^{(0)}, w_2^{(0)}, \cdots, w_d^{(0)})$の初期値を決める。(上付き文字は繰り返しの回数)
    3. $\mathbf{w}$がほとんど変化しなくなるまで↓を繰り返す。(要するに、$E$が最小になる$w$を一つずつ計算していく)
        1. $\frac{\partial E}{\partial w_1}(w_0, w_1, w_2^{(k)}, w_3^{(k)}, \cdots, w_d^{(k)}) = 0$になる$w_1$を求めて、それを$w_1^{(k+1)}$とする
        2. $\frac{\partial E}{\partial w_2}(w_0, w_1^{(k+1)}, w_2, w_3^{(k)}, \cdots, w_d^{(k)}) = 0$になる$w_2$を求めて、それを$w_2^{(k+1)}$とする
        3. $\cdots$
        4. $\frac{\partial E}{\partial w_d}(w_0, w_1^{(k+1)}, w_2^{(k+1)}, w_3^{(k+1)}, \cdots, w_d) = 0$になる$w_d$を求めて、それを$w_d^{(k+1)}$とする

### 座標降下法で使う数式

#### w0

$w_0$の計算方法

```math
E(\mathbf{w}) = \frac{1}{2} \sum_{i=1}^{n}{(y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij})^2} + \lambda \sum_{i=1}^{d} |w_i| 
```

より

```math
\begin{align}
\frac{\partial E}{\partial w_0} &= \frac{1}{2} \sum_{i=1}^{n} 2 (-1) (y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij}) \\
&= - \sum_{i=1}^{n} (y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij}) \\
&= - \sum_{i=1}^{n} (y_i - \sum_{j=1}^{d} w_j x_{ij}) + n w_0 \\
\end{align}
```

よって

```math
\begin{align}
\frac{\partial E}{\partial w_0} &= - \sum_{i=1}^{n} (y_i - \sum_{j=1}^{d} w_j x_{ij}) + n w_0 = 0 \\
w_0 &= \frac{1}{n} \sum_{i=1}^{n} (y_i - \sum_{j=1}^{d} w_j x_{ij})
\end{align}
```


#### w0以外

$w_k (k \neq 0)$の計算方法

```math
E(\mathbf{w}) = \frac{1}{2} \sum_{i=1}^{n}{(y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij})^2} + \lambda \sum_{i=1}^{d} |w_i| 
```

より

```math
\begin{align}
\frac{\partial E}{\partial w_k} &= \frac{1}{2} \sum_{i=1}^{n} 2 (-x_{ik}) (y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij}) + \lambda \frac{\partial |w_k|}{\partial w_k} \\
&= - \sum_{i=1}^{n} (y_i - w_0 - \sum_{j=1}^{d} w_j x_{ij}) x_{ik} + \lambda \frac{\partial |w_k|}{\partial w_k} \\
&= - \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij} - w_k x_{ik}) x_{ik} + \lambda \frac{\partial |w_k|}{\partial w_k} \\
&= - \sum_{i=1}^{n} \bigl( (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} - w_k {x_{ik}}^2 \bigr) + \lambda \frac{\partial |w_k|}{\partial w_k} \\
&= - \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} + \sum_{i=1}^{n} w_k {x_{ik}}^2 + \lambda \frac{\partial |w_k|}{\partial w_k} \\
\end{align}
```

よって$\frac{\partial E}{\partial w_k}=0$となる$w_k$は

```math
\begin{align}
\frac{\partial E}{\partial w_k} &= - \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} + \sum_{i=1}^{n} w_k {x_{ik}}^2 + \lambda \frac{\partial |w_k|}{\partial w_k} = 0 \\
\sum_{i=1}^{n} {x_{ik}}^2 w_k &=  \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} - \lambda \frac{\partial |w_k|}{\partial w_k} \\
w_k & = \frac {\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} - \lambda \frac{\partial |w_k|}{\partial w_k}} {\sum_{i=1}^{n} {x_{ik}}^2} \\
\end{align}
```


ここで

```math
\frac{\partial |w_k|}{\partial w_k}  = \left\{
\begin{array}{ll}
1 & (w_k > 0) \\
-1 & (w_k < 0)
\end{array}
\right.
```

を代入すると

```math
w_k = \left\{
\begin{array}{ll}
\dfrac {\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} - \lambda} {\sum_{i=1}^{n} {x_{ik}}^2} & (w_k > 0) \\
\dfrac {\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} + \lambda} {\sum_{i=1}^{n} {x_{ik}}^2} & (w_k < 0) \\
\end{array}
\right.
```

しかし、

```math
\begin{align}
\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} - \lambda < 0 \\
\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} + \lambda > 0 \\
\end{align}
```

すなわち

```math
- \lambda < \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} < \lambda \\
```

のとき、「$w_k > 0$となるべきだが、$w_k < 0$」「$w_k < 0$となるべきだが、$w_k > 0$」となってしまうので、このときは$w_k = 0$とする。よって

```math
w_k = \left\{
\begin{array}{ll}
\dfrac {\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} - \lambda} {\sum_{i=1}^{n} {x_{ik}}^2} & (\lambda \leq \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik}) \\
0 & (- \lambda < \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} < \lambda) \\
\dfrac {\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} + \lambda} {\sum_{i=1}^{n} {x_{ik}}^2} & (\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} \leq - \lambda) \\
\end{array}
\right.
```

となる。ここで、


```math
\mathbf{y} = 
\begin{pmatrix}
y_1 \\
\vdots \\
y_n
\end{pmatrix}
```

```math
\mathbf{I} = 
\begin{pmatrix}
1 \\
\vdots \\
1
\end{pmatrix}
```

```math
\mathbf{X} = 
\begin{pmatrix}
{x_1}^t \\
{x_2}^t \\
\vdots \\
{x_n}^t
\end{pmatrix}  = 
\begin{pmatrix}
x_{11} & \cdots & x_{1d} \\
\vdots &        & \vdots \\
x_{n1} & \cdots & x_{nd} \\
\end{pmatrix}
```

```math
\mathbf{x'} = 
\begin{pmatrix}
x_{1k} \\
\vdots \\
x_{nk}
\end{pmatrix}
```

```math
\mathbf{w_k'} = 
\begin{pmatrix}
w_1 \\
\vdots \\
0 \; (w_k = 0) \\
\vdots \\
w_d
\end{pmatrix}
```

とおくと、以下の行列計算でも表せる。

```math
w_k = \left\{
\begin{array}{ll}
\dfrac {(\mathbf{y} - \mathbf{I} w_0 - \mathbf{X} \mathbf{w_k'})^t \mathbf{x'} - \lambda} {\mathbf{x'}^t \mathbf{x'}} & (\lambda \leq \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik}) \\
0 & (- \lambda < \sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} < \lambda) \\
\dfrac {(\mathbf{y} - \mathbf{I} w_0 - \mathbf{X} \mathbf{w_k'})^t \mathbf{x'} + \lambda} {\mathbf{x'}^t \mathbf{x'}} & (\sum_{i=1}^{n} (y_i - w_0 - \sum_{j \neq k} w_j x_{ij}) x_{ik} \leq - \lambda) \\
\end{array}
\right.
```


# 図作ったときのコード置き場

- https://github.com/kzy611/qiita/blob/master/LinearRegression.ipynb

# 参考

- Coordinate descent(Wikipedia): https://en.wikipedia.org/wiki/Coordinate_descent